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ABSTRACT 


fr 


It  Is  commonly  known  that  In  non-honogeneous  media  (phase  velocity 
dependent  on  location)  refraction  of  acoustic  signals  occurs.  Solving 
the  wave  equation  with  variable  c  is  extremely  involved  and  the  cases 
where  solutions  can  be  found  do  not  give  very  rtuch  insight  into  the  physical 
meaning  of  the  problem. 

The  method  of  ray  tracing,  the  solution  of  the  eikonal  equation,  readily 
adapts  itself  to  non-homogeneous  media  and  describes  the  propagation  of 
wavefronts.  It  has  been  used  extensively  in  underwater  acoustics  but 
not  so  nuch  in  atmospheric  applications.  Some  reasons  for  the  limited 
use  of  ray  tracing  techniques  in  outdoor  sound  propagation  are  that 
1)  most  acoustic  work  in  recent  years  has  been  for  underwater  applications 
due  to  ITavy  sponsoring,  2)  and  also  that  very  few  simultaneous  measurements 
of  acoustical  and  meteorological  data  have  been  performed. 

Atmospheric  sound  ranging  techniques  have  in  the  past  neglected  verticie 
velocity  gradients.  Ray  tracing  is  a  useful  method  in  studying  propagation 
in  air  and  can  be  used  as  an  adjustment  to  sound  ranging  methods  to 
consider  atmospheric  variations. 

Presented  here  is  a  derivation  of  the  eikonal  equation  and  its  solution 
with  an  attempt  to  give  physical  reasons  for  this  approach.  A 
computer  nodel  of  the  technique  of  ray  tracing  for  atmospheric 
applications  (also  an  eigenray  nodel)  has  been  developed  and 
some  results  are  given  using  data  collected  in  field  measurements. 
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I*  Introduction 

The  concepts  presented  in  this  paper  are  by  no  means  new. 

It  is  hoped  that  the  approach  used  will  help  to  simplify  and 
clarify  a  study  that  has  been  complicated  by  mathematical  gymnastics.  The 
theory  presented  is  rigorous  but  the  steps  are  logical  and  it  is  not 
assumed  that  the  reader  is  already  familiar  with  ray  tracing. 

Ray  tracing  is  an  approach  that  was  developed  in  the  field  of 
optics.  Geometrical  optics*  as  it  is  called*  has  been  used  widely  in 
different  aspects  of  acoustics.  It  is  most  commonly  used  (in  acoustics) 
in  the  specialties  of  fluid  dynamics,  shock  theory  *  and 
non-linear  acoustics  and  called  the  method  of  characteristics. 

The  methodology  in  these  specialties  is  different  than  for 
sound  propagation  theory  but  the  approach  is  very  similar  and 
the  equations  take  the  same  form. 

Ray  tracing  techniques  have  been  used  for  many  years  in 

underwater  sound  propagation.  In  recent  years  many  acoustic 

approaches  have  been  used  In  meteorology.  Of  these  SODAR  (SOund 

Detection  And  Ranging)  has  been  used  to  determine  acoustic  rays 

and  from  resulting  data  to  approximate  temperature  profiles  under 

4 

inversion  conditions  (Increase  of  temperature  with  height). 

The  method  of  ray  tracing  has  been  promoted  in  the  field  of  outdoor 
sound  propagation  partially  due  to  new  Interest  in  noise  control. 

In  sound  ranging  applications  the  distance  to  the  sound  source  is 
different  than  simply  the  product  of  sound  speed  and  travel  time 
in  non-homogeneous  media.  Ray  tracing  is  seen  as  a  useful  method  in  the 
study  of  propagation  paths  in  non-homogeneous  media  where  refraction 
is  present. 

A  brief  pre»«n  stl  <  will  be  given  of  the  equations  leading  up 
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to  the  wave  equation*  The  eikonal  equation  will  be  derived  assuming 
a  series  solution  to  the  wave  equation  and  taking  the  first 
terms  of  the  expansion*  Solution  of  the  eikonal  equation 
will  be  first  done  in  the  homogeneous  (medium)  case  and  then 
in  the  non-hooogeneous  case.  Discussion  then  follows  concerning 
caustics  (high  concentration  of  energy)  and  shadow  zones  (zones 
of  silence)* 

A  general  discussion  of  the  computer  models  will  be 
presented  and  analysis  of  some  data  from  field  measurements 
will  be  analyzed  and  discussed*  For  more  information 
concerning  the  use  of  the  two  computer  programs  see 
Appendix  C.  Appendix  A  contains  the  listing  of  an  eigenray 
computer  program  which  solves  the  eikonal  equation  for 
rays  that  start  at  a  given  source  location  and  pass  through 
a  given  receiver  location.  Appendix  B  contains  a  ray  tracing 
routine  which  takes  source  location  and  starting  angles  either 
from  the  eigenray  program  or  from  some  other  source  and 
plots  the  resulting  rays. 

The  present  program  package  has  been  designed  to  analyze 
ray  paths  over  a  fiat  terrain  with  specified  vertical  temperature 
and  wind  profiles*  Attenuation  because  of  spherical  spreading 
and  atmospheric  absorption  has  been  included* 

Work  is  progressing  to  include  ground  effects  and 
variable  topography  in  the  package.  An  eigenray  routine 
designed  for  underwater  use  called  CONGRATS  (CONtinuous  Gradient 
RAy  Tracing  System)  is  being  revised  for  atmospheric  work. 

CONGRATS  fits  a  continuous  gradient  to  a  discrete  profile 
input.  The  final  routine  will  also  include  the  ability 
to  change  the  temperature  and  wind  profiles  in  a  path. 
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11.  Ray  Tracing  Theory 


A.  Derivation  of  the  wave  equation 


For  the  sake  of  completeness  the  place  to  begin  this  study 
is  with  the  basic  equations  leading  up  to  the  wave  equation* 

The  potential  velocity  4  is  defined  by 

V  *7  i  (1) 

The  approach  presented  here  is  constructed  around  the  potential 
velocity  but  it  Is  noted  that  it  can  be  developed  around  other 
quantities  such  as  velocity  or  pressure  equally  as  well* 

The  second  equation  needed  is  a  statement  of  Newton's 
first  law,  that  stress  is  equal  to  the  negative  of  momentum 
flux.  This  is  called  Euler's  equation  and  has  the  form 

VP  "  -p  3V/  3t 
o 

Substituting  equation  (1)  here  and  after  minor  manipulation, 
Euler's  equation  for  potential  velocity  becomes 


■  p  di  /  dt  +  constant 


(2) 


The  state  equation  is  a  statement  that  pressure  is  a 

function  of  density  If  expanded  in  a  series  around  the  ambient 

density  and  only  the  first  two  terms  are  retained  the  linearized 

state  equation  becomes 

p  ■  ?'( p  )  p 
o 

where  p  is  the  ambient  density* 
o 

◦f  the  propagation  speed.  So 

2 

P  *  c  p 


P'(p  )is  equal  to  the  square 
o 


(3) 


is  the  linearized  state  equation  to  be  used* 

The  final  equation  necessary  to  derive  the  wave  equation 
is  a  statement  of  conservation  of  mass  called  the  continuity 
equation.  It  states  that  the  net  flow  of  mass  into  (or  out  of) 
a  volume  Is  equal  to  the  net  change  of  mass  in  that  volume  and 
has  the  form 


4- 


-♦ 

7* v  -  -i/p  3 p/3 1 
o 

Substituting  from  equation  (3)  for  p  and  from  equation(l)  for  v, 
this  equation  becomes  ^ 

7  6  *  *  l/(poc  )  3p/3t 

And  substituting  for  p  from  Euler's  equation  (2)  the  wave 


equation  results 


2 

7  6 


2  2  2 
1/c  3  6/3 1 


(4) 


B.  Derivation  of  the  eikonal  equation 


The  eikonal  equation  is  a  transformation  of  the  wave  equation 
describing,  instead  of  the  wave  itself,  the  propagation  of  wave 
surfaces  or  wavel coats#  Rays  may  be  considered  as  packages  of 
acoustic  energy  travelling  normal  to  the  wavefronts#  Wavefronts 
are  the  loci  of  points  which  undergo  the  same  motion  at  a  given 
instant# 

Rays  in  this  theory  are  somewhat  equivalent  to  characteristics 
in  the  method  of  characteristics  used  in  both  non-linear  acoustics 
and  shock  theory*  The  difference  is  that  characteristics  take 
the  role  in  these  other  specialties  as  carriers  of  discontinuities. 
The  theories  are  very  closely  related.  In  section  C  a  solution  to 
the  eikonal  equation  is  developed  using  techniques  typical  of  the 
method  of  characteristics# 


To  derive  the  eikonal  equation  we  begin  by  defining  the  wavefront 
by  the  equation 

S(x, t)  -  0  (5) 

-♦AAA 

where  x*xi+xj+xk#  In  this  analysis  S  has  the  dimension 
12  3 

of  time  and  for  a  fixed  point  can  be  thought  of  as  the  difference 
between  the  time  that  has  past  and  the  time  necessary  for  the 
wave  defined  by  6  to  reach  that  point.  Therefore  for 

S(x,t)  <  0 


it  is  seen  that 


6(x, t)  -  0 
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Since  S  is  constant  in  reference  to  a  location  on  the  wave 
(e*g.  the  wavefront),  i  can  be  expanded  in  a  series  (Taylor) 
of  the  fora 


i 


S>0 

(6) 

S<0 


It  will  be  seen  later  that  0  (x)  represents  the  variation  in 

o 

magnitude  of  the  wave,  or  the  factor  of  spreading  loss. 

It  is  intended  chat  the  series  solution  will  be  substituted 

into  the  wave  equation  in  order  to  obtain  equations  for  S  and  0  • 

n 

The  series  of  equation  (6)  is  chosen  due  to  the  property  that 

when  derivatives  of  i  are  taken  the  derivatives  of 

H  -  Sn/n! 
n 

are  simply  H  ,  ie. 
n- 1 


H'(S)  -  H  (S) 
n  n-1 


Since  the  wave  equation  has  two  derivatives  the  form  of  H  for 

n 

negative  n  must  be  considered.  If  equation  (6)  is  rewritten  as 

00 

i  -  £  0  (x)  H  (S)  (7) 

n«0  n  n 


and  H  (S)  is  defined  as 
n 


H  (S) 
n 


S>G 

S<0 


it  is  noted  that  H  (S)  is  simply  the  Heaviside  function  and 
0 

the  H  (S)  are  simply  its  integrals.  Therefore,  using 
n 
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generalized  functions,  H  (S)  for  negative  n  can  be  defined  where 

n 


H  ^(3)  -  5 (S) 

and 

H  (S)  -  6 '(S) 

-2 

etc. 

Talcing  derivatives  of  equation  (7)  yields 


a* 

m 

°° 

I  e  (x)  h  (s) 

as 

ac 

n*0  n  n-l 

at 

CO 


U  > 
2 

7  6 


And  therefore  the  wave 


00 


L  (Ve  h  (s)  +  e  ii  (s)  7S) 

n*0  n  n  n  n- 1 

£  2 

■  X  V  9  H  (S)  +  2V0  •  (H  (S)  VS) 
n*0  n  n  n  n- 1 

2  2 

+  8  II  (S)  (VS)  +  0  H  (S)  V  S 
n  n-2  n  n-l 


equation  becomes 

.  2  1  /ds 

C6n  Hn-5(S)  ^7S)  "  /  — 

7\3t 


fr  2  1 

/  ^ 

r  2  I  a  s  1 

vs- 

0  +  2V0  •  VS  ► 

~2~~2 

n  n 

1  L  c  a  t  J 

V 

2 

+  v  e  ii  (s)  )  -  o 
n  n 


Grouping  like  terms  gives 


i 
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C .  Solution  to  the  eikonal  equation 

1.  General  Discussion 

To  get  a  general  feel  for  rays  and  what  the  eikonal  equation 
says,  a  perturbation  approach  is  in  order*  First  of  all 
che  unit  normal  to  the  surface  S  is  given  by 
a  -VS 

r  -  —  (10) 

;  vs  i 

Considering  an  initial  position  of  the  wavefront  depicted  by 

S(x  ,t  )  -  0  (11) 

0  0 

and  slightly  perturbing  all  of  the  variables  of  space  and  time 
the  surface  is  then  defined  by 

S(x  +  r  3,  t  +A  t)  -  0  (11a) 

0  0 

Then  a  derivative  may  be  approximated  by  a  finite  difference 
between  equations  (11a)  and  (11)*  The  result  is 

rVS  as  +3  S  At  *  0 

and  the  ray  velocity  (normal  to  the  wavefront)  is  then  given  by 

-as 

lim  £s  - 

A 

At*0  At  rVS 

A 

Substituting  for  r  from  equation  (10)  yields 

ds  -  36/ at  (12) 

dt  I  VS  I 

From  the  eikonal  equation  (8)  we  see  that 
ds/dt  *  +c 


"  4 
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Therefore  we  can  say  that  the  eikonal  equation  in  general 
says  that  a  wavefront  has  a  normal  velocity  or  the  ray  has 
a  velocity  of  magnitude  +c.  This  may  appear  to  be  a  trivial 
result  but  its  importance  is  that  this  does  correspond 
directly  to  the  solution  of  the  wave  equation.  It  means 
that  no  matter  what  changes  in  direction  a  wavefront 
may  undergo  it  will  propagate,  in  isotropic  media,  at  the 
characteristic  phase  velocity  of  the  medium  at  its  location. 

2.  Homogeneous,  isotropic  media 

By  homogeneous,  it  is  meant  that  the  propagation  velocity 
c  is  constant  with  respect  to  location  and  time  or  equivalently 
that  temperature*  is  constant  (isothermal  condition)  and  there 
is  no  wind.  Isotropic  conditions  imply  that  c  is  the  same 
regardless  of  the  direction  of  propagation.  This  is  one  of 
the  simplest  of  cases.  The  solution  shows  the  equivalence 
of  the  eikonal  equation  under  conditions  that  will  be 
demonstrated  later  to  the  wave  equation.  It  Is  common  to 
consider  homogeneous,  isotropic  media  when  solving  the  wave 
equation  but  the  power  of  the  ray  technique  is  seen  best  when 
these  conditions  are  relaxed. 

Specifying  the  wavefront  S  as 

A  ♦ 

S(x, C)  •  t  -  u(x)  -  0  (13) 

*♦ 

it  can  be  s.en  that  u(x)  locates  the  wavefront  at  x  for  varioua 
times*  Substituting  this  into  equations  (8)  and  (9)  yields 

-*  2  2 

(Mx)>  -  1/c  (I*) 

♦  2 

2?u(x)-?6  +  7  u  8  -  0 

0  0 


and 


(15) 
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Solving  equation  (14)  then  gives  a  solution  to  the  eikonal  equation* 
Consider  now  the  change  of  some  quantity  along  the  ray,  ie. 
the  first  derivative  in  terms  of  the  distance  s,  d/ds*  7h(x)  is 
normal  to  the  wavefront*  Equation  (14)  says  that  cVU(x)  is  unity 
and  therefore  this  represents  the  unit  normal  to  the  wavefront. 
Multiplying  this  quantity  by  the  change  along  x  once  again  yields 
d/d s  ie. 

d()  -  c7u'7()  (16) 

ds~~ 

This  equation  says  that  the  change  along  the  path  is  equal  to  the 
change  normal  to  the  wavefront.  The  so  called  characteristic  equations 
are  all  derived  directly  froa  equation  (16).  These  are  the  derivatives 
with  respect  to  s  of  x,  7u,  and  u.  Therefore 


dx  ■  c7u 


And  since  ?ii  is  constant  from  equation  (14) 


dVu  *  c7u  (7’(7u))  ■  0 


2  2 

du  ■  c(7u)  *  c/c 


Since  7u  is  normal  to  the  wavefront,  equation  (17)  shows  that  the 
rays  are  also  normal  which  is  how  we  initially  defined  rays.  It 
is  noted  that  in  anisotropic  media  the  rays  are  not  necessarily 
orthogonal  to  the  wavefront  (see  section  II-C-9)  .  Equation  (18) 
•ays  that  Vu  is  constant  along  the  ray.  It  is  concluded,  therefore, 
that  the  rays  are  straight  lines,  ie.  7u  is  constant  and  c  is 
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constant,  therefore  from  equation  (17),  the  path  along  the  ray  and 
the  vectorial  distance  vary  by  a  constant  and  the  ray9  oust  be 
straight  lines.  Equation  (19)  integrates  to 

u-  s/c  (20) 

which  means  due  to  equation  (13)  that  for  any  tine  t  >  0  that  the 
wave  surface  t  “  u  -  s/c  is  at  a  distance  ct  along  the  ray,  ie. 
s  »  ct.  These  equations  together  state  that  rays  can  be  constructed 
by  drawing  straight  lines  from  the  initial  wavefront.  Figures  1 
and  2  are  examples  of  this  construction  showing  a  spherical  source 
and  a  plane  source,  respectively. 


3.  Energy  conservation  and  attenuation  due  to  spreading 

It  is  noted  that  equation  (15)  can  be  rewritten  as 

2 

7-  (Vu  9  )  -  0  (21) 

0 

This  is  in  a  divergence  form  which  usually  indicates  the  conservation 

of  something.  It  is  common  to  think  of  this  as  an  equation 

showing  the  conservation  of  energy.  If  we  consider  a  flow  from  the 

wavefront  W  at  time  u  *  0  to  the  wavefront  W  at  time  u  »  t  as  shown 
1  2 
in  figure  3  and  Integrate  over  the  volume  defined  by  a  narrow  tube 

between  the  wavefronts  we  will  obtain  the  constant  energy  flux  law 

and  be  able  to  find  the  attenuation  due  to  the  rays  becoming  less 

dense  (le.  spreading  loss).  First  we  use  the  divergence  theorem 
6 

defined  by 
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Figure  1  -  Propagation  from  a  spherical  source  in  an  homogeneous 

isotropic  medium 


i 


-15 


fJF 


V 


F  dV 


A 

n  dS 


on  Che  volume  integral  of  equation  (21)  and  obtain 


2  a 

r ,  (Vu  9  )  •  n  -  0  (22) 

1)  0 
S 

A  + 

On  the  sides  of  the  narrow  tube  (see  figure  3)  n  and  Vu  are 
orthogonal  and  therefore 


A 

Vu*  n  *  0  on  the  sides. 

A 

On  tf  and  W  ,  n  and  Vu  are  in  the  same  and  opposite  direction 
2  I 

respectively,  therefore 


Vu* n  *■  |  Vuj  on 
»  -  |Vu|  on 


Also  from  Che  eikonal  equation  (14) 

|Vu|  -  1/c 

Therefore  equation  (22)  is  equivalent  to 


and 


2 

„  1  e 
Jf  Z  0 


dS 


0 


2 

,fi6  dS  -  0 
//  c  0 
W 

1 


and  since  c  is  constant  in  this  case  we  may  say 


2 

f f  0  dS  * 
0 
w 

1 


2 

ffe  as 
o 
w 

2 


(23) 


(24) 
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If  we  assume  that  we 
with  cross-sectional 
the  integrals  may  be 


are  integrating  over  a  small 
areas  aA^,  and  aA^  on  And 
approximated  by 


package  of  rays 
W  respectively, 


9  (x  )  aA  -  8  (x  )  aA 
0  1  l  0  2  2 


This  gives  in  the  limit  as  a  A  and  a  A  go  to  zero 

12 

-►  +  1/2  -1/2 
©  (x  )/9  (x  )  *  (dA  /dA  )  -  (dA  / dA  ) 

0  2  0  1  1  2  2  1 


(25) 


The  acoustic  ray  is  the  path  of  propagation  of  acoustic  energy. 
Equation  (25)  means  that  divergence  or  convergence  of  rays  indicates 
decreasing  or  increasing  energy  concentration,  respectively.  For 
example  in  plane  waves 

dA  /dA  -  1 
2  1 


ie.  the  cross-sectional  area  of  a  bundle  of  rays  stays  constant 
along  the  propagation  path  and  the  rays  are  parallel.  This  indicates 
that 


9  *  constant 

0 

or  that  there  is  no  spreading  loss.  For  cylindrical  and  spherical  rays 


dA  /dA  -  R 
2  1 

2 

dA  /dA  -  R 
2  1 


respectively. 


After  substituting  into  equation  (25)  we  have 
-1/2 


9  oc  ft 
0 


for  cylindrical  rays 


and 


-1 

9  a  a 
0 


for  spherical  rays. 


These  terms  are  consistant  with  spreading  losses  associated  with  wave 
phenomena.  In  a  non- homogeneous  medium  the  losses  will  be  similar. 

A  ratio  of  sound  speeds  at  one  wavefront  to  the  other  will 
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be  Included  (ie.  c  *  c(x,y,z))  but  we  may  assume  that  the 

ratio  will  be  near  I  and  therefore  can  use  the  previous  equations 

for  9  as  the  spreading  factor*  (See  section  III-A  for  a  more 
0 

precise  spreading  factor  dependent  on  range  rather  than  the 
distance  travelled). 

4.  Prediction  of  Caustics 

The  examples  presented  in  the  last  section  for  equation  (25) 
concerned  divergent  rays  and  therefore  showev  examples  of 
spreading  loss.  Another  effect  predicted  by  equation  (25) 
is  that  of  caustics.  Caustics  arise  when  an  initial 
wavefront  is  concave  away  from  the  direction  of  propagation 
causing  a  focusing  effect  as  in  figure  4.  The  cusp 
shaped  envelope  is  called  a  caustic.  The  region  inside 
the  envelope  is  triply  covered  by  rays  and  energy  is 
concentrated.  On  the  caustic  neighboring  rays  touch  each 
other  and  therefore  the  bundle  of  rays  described  in 
the  last  section  has  a  cross-sectional  area  of  zero  ie. 

dA  /dA  -  0 
2  1 

which  predicts  from  equation  (25)  that 

e  -  00 

0 

Caustics  or  points  In  space  where  there  is  infinite^ 
acoustic  energy  are  also  predicted  by  the  wave  equation. 

The  question  here  is  whether  the  linearized  wave  equation  (4) 
applies  in  this  case.  It  should  be  recognized  that  at 
caustics  there  is  high  acoustic  energy  concentration  but 
because  of  non-linear  effects  it  Is  not  infinite. 

Caustics  will  also  be  evident  In  non-homogeneous  and 
anisotropic  media  but  are  not  as  easily  described  as  the 
cusp  shaped  envelope  which  arises  in  homogeneous,  isotropic 
media. 
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5,  Non-homogenous ,  isotropic,  stratified  media 


By  n on-homogeneous  it  is  meant  that  the  phase  velocity 
depends  on  location  ie. 

c  -  c(x,yfz) 

In  most  cases  c  is  considered  as  a  function  of  height  only, 
but  for  comparison  the  more  general  equations  are  presented 
and  then  the  equations  for  the  simpler  stratified  case* 

The  only  difference  in  the  eikonal  equation  (14)  is 
that  c  is  no  longer  constant.  Using  eqvtion  (16)  will  still 
give  the  proper  characteristic  equations  for  dx/ds,  d  yu  /ds 
and  du/ds  as 

dx  -  cyu  (26) 

ds 

2 

dyu  *  cyu(y*  yu)  *  1  cy  *( yu) 
ds  T 


and  because  of  the  eikonal  equation  (14) 

-2 

dyu  *  L_  cy*  (c  ) 
ds  2 


-c 

~T 


yc 


c 


therefore 


dyu  « 
ds 


and 


2 

du  *  c(yu) 
ds 


c  *  1 
“T  ~ 
c  c 


(27) 


(28) 


Since  yu  is  normal  to  the  wavefront  equation  (26)  like  equation  (17) 
says  chat  the  rays  are  also  orthogonal  to  the  wavefront. 
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However,  equation  (27)  is  different  than  equation  (18)  and 
says  that  cyu  is  no  longer  constant  on  a  ray  path  and  the 
combination  of  these  two  equations  says  that  the  rays  bend 
around  in  response  to  the  gradient  of  the  phase  velocity  (tc). 

The  negative  sign  in  equation  (27)  indicates  that  the  rays  bend 
toward  a  region  of  lower  velocity.  Solving  equations  (26)  and  (27) 
sirault aneously ,  gives  the  rays  and  equation  (28)  gives  the  travel 
time  by 


u  *  r  ds 

J  cO 


along  ray 


(x) 


(29) 


In  a  stratified  medium  these  equations  simplify  considerably 
so  that  they  are  more  easily  solved*  In  this  case  the  phase 
velocity  depends  only  on  height  ie. 

c  -  c(z) 

The  characteristic  equations  become 


dx  *  c7  u 

dz 

x.y 

ds 

ds 

d  v  u 

*.y  ”  o 

ds 


(30) 


(31) 


du  *  1^ 
ds  c 


(32) 


where  x  and  7  are  the  horizontal  components  and  gradient 
x.y 

respectively,  and  z  and  7  are  the  verticle  component  and 

z 

gradient.  From  equation  (31)  7  u  is  constant  and  from 

x.y 
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equation  (30)  we  see  that  dx/ds  gives  the  angle  from  the 
horizontal  that  the  ray  makes  as 


dx  »  cos  0  •  c(z)  7  u 

. 

ds 


If  a  subscript  zero  refers  to  an  initial  point  we  have 

cos  0 


7  u 


(33) 


and  since  7  u  is  constant  we  can  write  the  equation 

x,y 

cos  0  *  c(z) 

cos  0  c 
0  0 


(34) 


which  Is  Snell's  law  in  optics.  From  the  eikonal  equation 
we  know  that 

2  2  2  2 
(7  u)  +  (7  u)  -  (7u)  -  1/c 

x,y  z 


and  substituting  from  equation  (33)  gives 


I- 


1 


cos  0 


l  C  (z) 


ol 


1/2 


C  J 
0 


(35) 


to  solve  for  rays  the  ray  equations  (30)  may  be  combined  into 
dx  •  7  u 

nr  -HV 


Substitutlng  equations  (33)  and  (35)  and  Integrating  yields  the 
equation 


*  ■*  /  c(z)  cos  8  /c 

x  -  x  "  /  0  0 

0  /  - Z - 2 - r 

J  (1-c  (z)cos  0  /c  ) 
z  0  0 


m 


dz 


(36) 
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which  describes  a  ray  with  initial  angle  0^  at  a  point  (x^,y^,z^)« 

From  equations  (12)  and  (10)  the  ray  travel  time  is  given  by 

z 


u  -  f  ds  -  C  dz 
J  c  /  c  V  u 


or 


u  *  r _ 

J  c(z) 


dz 


7 - 2 - 7rm 

(1  -  c  (z)  cos  0  /c  ) 

0  0 


(17) 


These  last  equations  (36)  and  (37)  are  the  basis  of  the  model 
presented  in  this  paper.  In  the  next  section  the  question  of  when 
the  eikonal  equation  is  valid  is  considered  followed  by  a  section 
which  discusses  two  particular  phase  velocity  distributions. 


b.  Conditions  of  validity  of  the  eikonal  equation 

It  is  emphasized  that  the  eikonal  equation  is  only  an  approximation 
to  the  linearized  wave  equation.  The  word  linearized  is  stressed  so 
that  one  is  aware  that  the  wave  equation  itself  is  not  always  valid 
and  certainly  an  approximation  to  it  would  not  be  valid  under 
non-linear  conditions.  One  of  these  conditions,  that  resulting 
in  caustics,  has  already  been  noted  in  section  II-C-4. 
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In  chis  section  an  harmonic  solution  to  the  wave  equation  (4)  is 
considered  and  substituted,  resulting  in  the  eikonal  equation  with 
some  descrepancy*  Making  this  error  small  is  the  condition  sought, 
so  that  the  eikonal  equation  will  be  a  good  approximation  to  the  wave 
equation* 

First  the  assumption  is  made  that  the  solution  is  time  harmonic 
only  if  the  wave  has  reached  the  spatial  coordinate  specified*  The 
wavefront  S  as  defined  in  equation  (13)  is  an  appropriate  time 
frame  to  consider*  It  is  also  assumed  that  the  amplitude  of  the 
wave  may  vary  In  space  due  to  variations  in  the  medium*  The 
solution  is  then  of  the  form 

4  4 

i  »  A(x)  exp(>S(x,t)) 

*4  *4> 

*  A(x)  exp(jw(t  -  u(x))) 

Substituting  into  the  wave  equation  (4)  the  resulting  equation  is 

2  2  2  2  2  2 
V  A  -  v  A(Vu)  -  J(2wVa*Vu  +  o>AV  u)  *  -w  A/c 

Separating  the  real  and  Imaginary  parts  yields 

2  2 

-  1  V  A  +  (Vu)  -1^-0  (38) 

ai  A  c 

and 

2 

V  u  +  £  VA-Vu  -  0  (39) 

A 

For  u  to  be  a  solution  to  the  eikonal  equation  the  first  part 
of  equation  (38)  must  be  zero*  This  will  be  so  if  the  amplitude 
of  oscillation  A  is  constant  or  linear  in  which  case  the  second 
spatial  derivative  of  A  would  be  zero;  or  if  the  frequency 
is  infinite*  In  general  neither  of  these  assumptions  can  be  made* 

The  previous  condition  may  be  relaxed  by  making  the  first  term 
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in  equation  (33)  much  smaller  than  the  second  ie. 


2  2 
7  A  «  (7u) 

<i>  A 


Using  the  elkonal  equation  (14)  this  becomes 

2  2 


2  2 
c  7  A 
7 


V  2tt  /  A 


«  1 


(40) 


u  A  V  2*  /  A 

for  convenience  the  gradient  of  a  function  will  be  defined  over 
the  distance  of  one  wavelength  so  that 


7F  -  AF/\ 


(41) 


This  will  transform  equation  (40)  into 

x  A  7 A  «  1  (42) 

A 

If  this  condition  is  met,  u  is  a  solution  to  the  elkonal  equation* 
For  u  to  also  be  a  good  approximation  to  the  wave  equation  or 
rather  for  the  ray  solution  to  be  a  good  approximation  to  the 
wave  solution,  equation  (39)  must  also  be  satisfied  or 

2 

7  u  -  -2  7A  .  7u  -  -2  7A 

A  Ac 


and  from  equation  (42)  this  gives 

2 

-  X  A  c7  u  «  1  (43) 


Taking  the  gradient  of  the  elkonal  equation  (14)  we  have 

2  -3 

2(7u)  7  u  -  -2  c  7c 


or  using  the  square  root  of  the  elkonal  equation 

2 

7  u  ■  -  7c 
T 
c 
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Substltutlng  Into  equation  (43)  we  have 

,  *A  «  1 

c 

or  from  equation  (41 )  and  knowing  that  changes  in  c  are  small 

A  Vc  «  Vc  (44) 

This  is  the  condition  sought.  It  states  that  a  solution  to  the 
elkonal  equation  will  be  a  good  approximation  to  a  solution  of 
the  wave  equation  if  the  change  in  the  gradient  of  the  phase 
velocity  over  a  wavelength  is  small  compared  to  the  gradient 
itself.  Therefore  ray  solutions  are  valid  if  there  are  no 
large  changes  or  discontinuities  in  the  phase  velocity  profile. 

7*  Formation  of  a  shadow  zone  in  a  stratified  medium 

It  was  shown  in  section  II-C-5  that  rays  bend  toward  areas 

of  lower  phase  velocity.  From  this  simple  concept  it  can  be 

surmised  that  if  there  exists  a  maximum  phase  velocity  at  some 

height  z  above  a  source  that  a  shadow  zone  will  be  formed, 
m 

A  shadow  zone  is  an  area  where  no  acoustic  energy  penetrates. 

*More  specifically,  rays  near  the  height  z  will  bend  either 

m 

upward  or  downward  away  from  that  height.  This  is  illustrated 

in  figure  5.  The  limiting  ray  which  defines  the  boundary  of 

the  shadow  zone  is  that  ray  which  becomes  horizontal  at  height 

z  •  This  ray  is  often  called  a  split  -  beam  ray  since  it  may 
m 

be  bent  upward  or  downward  and  theoretically  is  handled  by 
considering  that  it  goes  both  ways. 

From  Snell's  law,  equation  (34)  we  have  that 
-1  -1 

8  *  cos  (c  cos  8  /c  )  ■  cos  (c  cos  8  /c  )  (45) 

0  0  mo 

So  that  as  c  increases  le  as  the  ray  nears  height  z  ,  8  will 

m 

decrease.  So  that  this  equation  is  defined  for  all  values  of  8 
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and  c,  c  can  be  chosen  Co  be  Che  aaxiaua  phase  velociCy  and  9 
a  a 

Che  angle  at  Che  height  associaced  vich  c  • 

a 

In  general,  if  Che  inicial  angle  of  a  ray  is  zero,  equation 
(45)  says  chat  Che  ray  will  not  be  horizontal  again  until  it 

reaches  a  level  with  c  •  c  .  In  figure  5,  if  Che  inicial  angle 

0 

-1 

9  >  cos  c  /c  (46) 

0  0  a 


then  9  will  be  greater  chan  zero  and  Che  rays  will  penetrate  Che 

a 

level  of  aaxiaun  sound  speed  and  continue  upward.  However,  if 

-1 

9  <  cos  c  /c  (47) 

0  0  a 


Che  cos  9  increases  Co  1  and  9  decreases  to  zero  at  Che  height  of 

z  defined  by  Che  equation 

c(z)  -  c  /cos  9  (48) 

0  0 


At  this  point  Che  ray  bends  downward. 

The  critical  ray  is  the  split  -  beaa  ray.  This  results  when 
9-0 


which  occurs  when 

-1 

9  -  cos  (c  /c  )  (49) 

0  0  a 


A  critical  distance  along  Che  ground  aay  be  defined  where  Che 
split  -  beaa  ray  intersects  Che  ground.  It  is  instructive  to  use 
a  case  as  in  figure  5  where  Che  velocity  gradient  is  linear  up  to  a 
aaxiaun  velocity.  Equation  (3b)  gives  us  Che  distance  travelled 
in  integral  fora  Placing  the  sound  source  on  the  ground  sets  z^  «  0 
and  we  can  define  x  -  0.  The  phase  velocity  is  then 


c  -  c  +  c  az 

0  0 


ie.  there  is  a  linear  velocity  gradient*  Equation  (36)  then 
becomes 


/  A ( 1  +  a z)  dz 

I 

J  (1  -  A  (A  +  az)  ) 


Where  A  ■  cos  0  •  Setting  r  *  1  +  az  and  dr  -  adz  gives 


l+az 

/  A  r  dr 

/  - TTH 

J  a(l  -  A  r  ) 


2  2  1/2 

-  1  (-<1  -  A  r  ) 


After  slight  manipulation  this  reduces  to 


2  1/212  f  12 

-  (1  -  A  )  j  +  jz  +  1}  -  1 


-J  "7-2 
a  a  A 


which  is  an  equation  for  a  circle*  Substituting  for  A  gives 

2  2 

f  x  -  tan  0',  +  f  z  +  1  ^  ■  1 

Oil  _ 


1  I 


a  cos  0 


which  is  a  circle  centered  at 


with  radius 


(tan  0^/a,  -1/a) 


R  *  1/a  cos  0 
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Thls  means  that  for  a  linear  velocity  gradient  rays  will  follow 
the  arc  of  a  circle* 

For  the  distance  travelled  along  the  ground,  z  *  0  and 


-  1 

"2 

a 


2 

tan  6 

0 


Therefore  the  ray  travels  horizontally 


x 


2  tan  8 

0 


a 


(50) 


before  reaching  the  ground  again*  Using  equation  (49)  for  O^defines 
this  distance  for  the  critical  ray  by 


x 

c 


-1 

2  tan  (cos  (c  /c  )) 
_ Q  m 

a 


(50a) 


where  a  is  the  slope  of  the  velocity  profile* 

Actually  rays  from  the  source  may  penetrate  the  shadow  zone  by 

multiple  reflections  off  the  ground*  If  in  addition  there  exists 

a  local  maximum  characteristic  velocity  above  the  maximum  c  , 

m 

rays  may  again  be  bent  downward  into  the  shadow  zone*  For  a  more 

precise  treatment  one  must  Include  the  effects  of  diffraction 

which  are  not  readily  defined  using  rays*  However,  the  existence 

of  shadow  zones  has  been  experimentally  observed  as  low  Intensity 
3,8 

zones* 
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3*  Waveguides 

Waveguides  result  from  the  existence  of  a  raised  minimum 
velocity.  This  is  illustrated  in  figure  6.  The  derivation  of 
this  result  is  similar  to  that  for  the  shadow  zone.  If  the 
initial  angle  is  specified  by  equation  (46)  the  ray  will 
penetrate  into  the  region  of  higher  velocity.  However  if 
equation  (47)  describes  the  initial  angle  the  cos  0  will 
increase  to  1  and  0  will  decrease  to  0  and  the  ray  will 
bend  downward.  At  this  point  the  ray  crosses  the  minimum 
value  of  c  again  and  will  repeat  the  pattern  symmetrically 
about  the  height  of  the  minimum. 

Waveguides  are  important  in  a  discussion  of  ray  theory 
since  it  allows  a  ray  to  propagate  for  a  long  distance 
without  reflections  and  probable  losses  from  ground 
interactions. 

9.  Anisotropic,  homogeneous  media 

In  anisotropic  media  the  phase  velocity  is  dependent 
on  orientation.  A  simple  example  is  when  sound  propagates 
in  a  wind.  The  sound  speed  will  be  greater  in  the  direction 
of  the  wind  than  orthogonal  to  it.  The  eikonal  equation  (14) 
still  remains  in  the  same  form,  with  the  phase  velocity 
now  a  vector,  ie. 

2  2 

(7u(x))  -  l/(c)  (51) 

-► 

In  this  case  only  homogeneous  media  are  being  considered  so  c  is 
constant.  The  characteristic  equations  arise  from  equation  (16) 


i 
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in  the  following  font, 
dx 


c(Vu-I) 


ds 

dVu 

-  0 


ds 

du  dx  *► 

*  7u»  *  7u*c(7u*I) 


(52) 

(53) 

(54) 


ds  ds 

Equation  (53)  says  that  the  gradient  of  u  is  constant  making  the 
rays  straight  lines  as  would  be  assumed  in  a  homogeneous  medium. 
However,  the  ray  direction  specified  by  equation  (52)  will  be 
parallel  to  the  wavefront  normal  if  and  only  if 

c  «  Vu 

This  will  be  true  if  and  only  if 


+  2  2  2  2 
(c*Vu)  -  F(p  +  p  +  p  ) 
12  3 


(55) 


where  the  p  's  are  the  components  of  7u  and  F  specifies  some  function. 
Equation  (55)  says  that  the  rays  are  normal  to  the  wavefront  only 
in  isotropic  media. 

Integrating  equation  (54)  and  substituting  u  *  t  locates  the 

wavefront  at  successive  times  ie. 

-►  ■+ 
u  *  s7u*c(Vu-I) 

or 

s  ■  (56) 

-t  ” 

7uc(Vu  I) 

The  vectorial  distance  is  specif  ^d  by  integrating  equation  (52)  as 

x  ■  sc(7u*I)  (57) 

From  equations  (51)  and  (57)  the  unit  vector  in  the  direction  of 
the  ray  is 


c(7ul) 
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Equation  (12)  gives  the  velocity  normal  to  the  wavefront  as 

ds  1  1 

_  *  _ »  _ -  c 

dt  (7s)  IVul 

Therefore  the  unit  vector  normal  to  the  wavefront  is 

cVu 

The  angle  m  between  the  normal  and  the  ray  path  can  then  be 
given  by 

cos  m  *  c(  7u- 1)  *  (c7u)  (58) 

Using  just  the  individual  components  of  the  vectors  in  this 
expression  will  give  the  angles  in  each  plane*  Using 
expressions  (58)  and  (56)  the  distance  along  the  ray  may 
be  specified  by 

ct 

s  -  _  (59) 

cos  m 

This  equation  says  that  the  wavefront  moves  along  the  ray  with 
speed  c / cos  m  which  is  greater  than  c. 


-I***' ’ 
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Iil-  computer  nodel  Description 
A.  ?asic  development  of  equations 

In  the  previous  sections  ray  -  theory  has  been 
developed  and  discussed.  This  section  is  devoted  toward 
techniques  used  in  the  developnent  of  computer  programs 
from  the  equations  derived.  Ihere  are  two  types  of 
programs  to  be  described.  These  are  1}  graphic  ray 
tracing  programs  and  2)  eigenray  programs.  In  both 
types  of  programs  the  sound  velocity  profile  must  be 
spec  if led • 

Since  the  sound  velocity  as  a  function  of  height  is  not 
easily  measured  other  related  units  must  be  measured  .  The 
sound  velocity  is  directly  proportional  to  the  square 
root  of  absolute  temperature  as  given  by 

1/2 

c  -  20. C5  (T) 

where  c  is  in  meters  per  second  and  T  is  in  degrees  Kelvin 

o 

(*  degrees  Celsius  +  273.2).  Since  this  refers  to  propagation 
relative  to  the  medium  we  must  include  the  wind  velocity  in 
this  forrulation  so  that  the  equation  specifies  propagation 
relative  to  the  ground  ie. 

1/2 

c  -  20.05  (T)  +  -7V  (60) 

Trie  factors  T  and  TV  can  be  measured  using  thermistors  and 
aneroneters,  respectively  and  the  vectorial  direction  of  the  wind 
using  a  hi  -  vane.  Therefore  the  phase  velocity  as  a  function 
of  height  may  be  specified. 

In  the  development  of  the  characteristic  equations  it 
was  necessary  to  use  the  vertical  phase  velocity  gradient  given 
by  dc/dz.  In  modeling  techniques  it  is  usual  to  use  a  linear 


35- 


difference  approximation  to  derivatives,  therefore 

*  Cl^“l  °i  (61) 

z  -  z 
i+1  i 

where  is  the  gradient.  Using  many  segments  for  the  gradient 
will  approximate  a  smooth  curve  fairly  well  and  therefore  other 
difference  forms  (e.g.  logarythmic)  are  not  used.  It  Is 
intended,  however,  that  for  a  small  number  of  values  of 
T  and  WV,  to  include  equations  from  meteorological  theory 
to  interpolate  other  values.  The  present  model  does  not  Include 
these  interpolation  methods. 

The  assumption  for  the  model  is  that  instead  of  a  simple 
stratified  medium,  the  medium  is  divided  into  layers  and 
each  layer  has  a  linear  gradient.  We  can  therefore  use 
the  equations  developed  earlier  to  derive  equations  for  each 
layer  and  follow  Individual  rays  from  layer  to  layer. 

Three  cases  must  be  considered:  i)  the  Isovelocity  case, 

2)  variable  velocity  when  the  ray  penetrates  the  layer  and 

3)  variable  velocity  when  the  ray  is  refracted  back  towards 
its  entry  level.  The  isovelocity  case  Is  really  simply  the 
homogeneous  case  discussed  in  section  II-C-2.  In  this  case 

g  *  0  and  the  rays  are  straight  lines.  If  D  is  the  thickners 

of  the  layer  and  9  is  the  angle  of  the  ray  upon  entering  the 

layer  the  change  in  the  x  distance  will  be  defined  by 

DX  -  D  cot  9  (62) 

i 

From  equation  (20)  the  travel  time  is  given  by 

2  2  1/2 
DS  (DX  *  D  ) 

DT  -  -  _  (63) 


c 


c 
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DT 


1 

■  In 
g 


1  +  sin  0(z) 
cos  9( z) 


1  (1  +  sin  0  )(i  -  sin  0  ) 

DT  -  _  in _ i _ i+1  (65) 

2g  (1  +  sin  0  )(i  -  sin  0  ) 

i+1  i 

The  third  and  final  case  is  when  a  ray  is  bent  around  and  returns 
in  the  direction  it  entered  the  layer.  First,  it  is  noted  from 
equation  (48)  that  if  the  ray  becomes  horizontal  at  a  point  where 
the  phase  velocity  is  given,  the  highest  value  of  z  is  defined  by 


Second,  it  was  shorn  in  section  II-C-7  that  rays  travel  in  a 
circular  path.  Also,  the  ray  may  turn  before  reaching  the  edge 
of  a  layer.  Therefore,  since  there  is  circular  motion,  the 
height  attained  in  a  layer  is  given  by 

1 

OZ  -  (l  -  cos  0  )  (66) 

—  i 

gfc 


where  0  and  z  are  measured  at  the  entrance  to  the  layer.  If 
i  i 

this  difference  is  greater  than  the  thickness  of  the  layer,  the 
ray  will  not  be  bent  around  in  that  layer.  If  DZ  in  this  equation 
is  less  than  D  then  OX  is  defined  by  equation  (50) 

2  tan  0 

OX  -  _ i  (50) 


a 
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where  a  is  Che  slope  of  Che  gradient  given  by 

g  *  c  a 
i  1 

The  time  DT  can  Chen  be  found  using  equaCion  (37)  with  liraics  of 

incegracion  z  and  z  +  D2  and  doubling  Che  resulc  since  Che  ray 
i  i 

muse  reCurn  Co  iCs  encry  height.  Therefore, 

1  I  -  sin  9 

DT  -  _  In  _ i  (67) 

g  1  +  sin  9 

i 

Equacions  (62),  (63),  (64),  (65),  (66),  (50)  and  (67)  form 
Che  basis  of  che  computer  models*  The  total  horizontal  distance 
and  time  Che  ray  undergoes,  x  and  C,  are  found  by  adding  all  Che 
DX's  and  DT's,  respectively*  The  actual  distance  the  ray  travels,  s, 
is  given  by  the  sum  of  DS's  where 

2  2  1/2 

DS  -  (D  +  DX  )  (68) 

for  the  homogeneous  case,  or  because  the  radius  of  curvature  is 
defined  by 

ds 

_  -  R 

de 

and  R  was  given  in  section  1I-C-7  as 

1  cos  9 


a  cos  9  c 

i  i 

therefore  for  the  non-horaogeneous  case 

cos  9 

DS  -  g  i  (9  -  9  )  (68a) 

i  -  i+1  i 

c 

i 

In  the  case  of  atmospheric  sound  propagation  there  are  only  reflections 
from  the  ground*  Ground  reflections  are  specular  and  handled  by 
taking  the  negative  of  the  angle  of  incidence* 
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The  first  type  of  program,  graphic  ray  tracing,  is  constructed 
from  these  equations  and  includes  the  reflections.  The  remainder 
of  rhls  program  consists  of  graphics  techniques. 

Input  to  the  ray  tracing  program  Includes  the  temperature  and 
wind  profiles  and  the  location  and  angle  of  the  sound  source.  For 
rays  travelling  upward  it  is  also  necessary  to  include  a  maximum 
height  that  is  to  be  considered.  This  height  can  sometimes  be 
conveniently  chosen  just  above  a  raised  inversion,  (velocity  is 
greater  at  a  greater  height).  Appendix  B  contains  a  graphics  ray 
tracing  model* 

Eigenray  routines  find  rays  that  travel  from  a  source  location 
to  a  specified  receiver  location.  This  is  accomplished  by  searching 
a  range  of  angles  and  using  a  bisection  method  to  zero  in  on  the 
angle  at  the  source.  The  program  follows  many  rays  by  the  method 
used  for  ray  tracing  and  internally  varies  only  the  source  angle 
until  a  solution  Is  found.  Once  this  is  completed  the  sound  field  at 
the  receiver  may  be  ascertained. 

In  the  prediction  of  the  sound  field  one  must  Include  the  effects 

of  absorption  and  spreading  losses.  To  obtain  the  intensity  spreading 

loss  a  solid  angle  ft  is  defined  with  symmetry  around  the  z-axis  so  that 

d  Q  m  2 n  cos  0  d0 
0  0 

where  the  angles  are  specified  in  figure  7.  The  unit  of  intensity 
will  be  defined  by  the  ratio  of  d U  to  the  area  dA  swept  out  by  the 
wave  surface.  From  figure  7  this  is 

dft  2*  cos  8  d0 
i  -  0  0 


dA  2 tr  x  sin  0  dx 
h 

cos  0  d0 

-  0  0  (69) 


x  sin  0  dx 
h 


*  trifCT - 
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The  horizontal  range  x  is  a  function  of  the  height  h  and  initial 

angle  6  therefore,  the  horizontal  unit  of  range  is 
0 

ax 

dx  -  de 
—  0 
36 

0 

Substituting  this  into  equation  (69)  and  taking  the  reciprocal  of 
the  resulting  function  will  yield  the  loss*  On  a  log  scale  this  is 

x  sin6  ax/36 

L  »  10  log  _ h _ 0  (70) 

cos  6 

0 


Equation  (64)  is  used  to  find  an  expression  for  ax/36^.  It  was  said 
that  x  is  the  sum  of  the  DX's,  therefore 


3  x 
36 


c  sin  6 
0 

- r~ 

cos  6 


2  23 


sin  6  -  sin  6 

i  i+1 


i-0 


c 

+  JO _ 

cos  6 

0 


36 

-  cos  6  1+1 

i+1  36 

0 


If  we  differentiate  Snell's  law,  equation  (34)  we  have 

3  6  c  sin  6  sin  6  cos  6 

_ i  -  J.  _ 0  -  0  i 

3 6  c  sin  6  cos  6  sin  6 

0  0  i  0  i 
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Substitution  into  equation  (71)  and  after  slight  manipulation  we  have 

n 

3  x  c  sin  8  l/l  1 

-  0  0  VT 

—  — 2 —  >  - 
d8  cos  8  \sin  8  sin  8 

0  0  i  \  i  i+1 

i-0 


sin  8 


cos  8 


■  DX. 

o  1 

^Hsin  8 


(72) 


10 


i-0 


sin  8 
i  i+1 


using  equation  (64).  Therefore  the  intensity  spreading  loss  is 

n 

x  sin  8  sin  8  DX 

L  -  10  log  _ ^n+1  0  ^ _ i  (73) 


cos  8 


i-0 


sin  8  sin  8 

i  i+1 


To  this  value  the  ground  absorption  and  atmospheric  absorption 
must  be  added* 

Presently  ground  losses  are  handled  simply*  The  number  of 

ground  reflections  n  is  counted  and  multiplied  by  a  loss  coefficient, 
b 

L  ,  provided  by  the  user*  It  is  intended  to  revise  this  by  using 
b 

a  closed  form  where  the  impedance  of  the  ground  will  be  specified 

and  phase  information  will  be  retained* 

The  atmospheric  absorption  coefficient  is  calculated  using 

the  American  National  Standard*  The  necessary  equations  are 

included  here  for  easy  reference.  The  absorption  coeffeicient  is 

2  -11  1/2 
Alpha  -  f  (1.84  X  10  (T/T  ) 

0 

-5/2  -2 

+  (T/T  )  (1.278  X  10  (exp(-2239.2/T)) 

0 

2  -1 

/(f  +  (f  /f  ))  +  1.068  X  10  (exp(-3352/T)) 

r,0  r,0 

/<f  +  (f  /f  ))))  (74) 

r,N  rfN 

in  Nepers  per  meter.  In  this  equation  T  is  the  temperature  in 
degrees  Kelvin  and  T^  is  the  ambient  temperature  equal  to  293,15  K; 


f  is  the  frequency  of  the  source  in  Hertz  and  f  and  f  are 

r,0  r,N 

the  relaxation  frequencies  in  Hertz,  for  oxygen  and  nitrogen 
respectively  ,  and  are  given  by 

4 

f  -  (24  +  4.41  X  10  h  X  ((0.05  +  h)/(0.391  +  h))) 

r  >0 

and  (75) 

-1/2  -1/3 

f  -  (T/T  )  (9  +  350h  exp  (-6.142((T/T  )  -  1))) 

r,N  0  0 

In  all  of  these  equations  the  pressure  Is  considered  equal  to  the 

ambient  pressure  and  so  doesn't  enter  Into  the  calculations.  For 

the  model  the  average  value  of  temperature  is  used  for  T. 

The  variable  h  is  the  per  cent  humidity  and  can  be  calculated  as 

h  -  h  (p  /p  )  (76) 

r  sat  so 

where  h  is  the  relative  humidity  and  the  ratio  of  saturation 
r 

pressure  to  ambient  pressure  can  be  calculated  from 

log  (p  /p  )  -  10.79586  (1  -  (T  /T)) 

10  sat  so  01 

-  5.02808  log  (T/T  ) 

10  01 

-4  -8.29692( (T/T  )-l) 

+  1.50474  X  10  X  (1  -  10  01 

-3  4.76955(1-(T  /T)) 

+  0.42873  X  10  (-1  +  10  01  ) 

-  2.2195983  (77) 


where  T  *  273.16  is  the  triple  point  isotherm  temperature. 

The  total  loss  is  then  given  by 

n 

x  sin  6  sin  8  DX 

TL  -  10  log  _ 0  n+1  ^  1 _  (78) 


cos  8 


sin  8  sin  8 

i  i+1 


+  Alpha  (x)  +  n  L 

b 

We  now  have  the  basis  for  an  eigenray  routine.  Appendix  A 
contains  such  a  model.  To  graph  the  eigenrays,  the  output  of  the 
program  in  Appendix  A  is  input  into  the  program  in  Appendix  B. 
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1#  Eigenray  routine  Improvements 

Since  the  present  model  is  for  a  horizontally  homogeneous 
medium  it  can  be  surmized  that  after  ground  reflections  and 
rays  reach  the  initial  height  and  angle  the  rays  will  follow 
the  same  pattern-  Advantage  is  taken  of  this  cyclic  nature  to 
speed  up  the  calculation  process.  It  is  necessary  to  calculate 
only  one  cycle  and  compare  the  horizontal  length  of  the  cycle  to 
the  range. 

Two  types  of  intersection  with  the  receiver  are  possible 
within  one  cycle;  1)  as  the  ray  is  upward  bound  and  2)  as  the  ray 
goes  downward*  A  range  of  initial  angles  is  swept  through  and 
rays  coming  near  the  receiver  location  are  stored. 

The  rays  then  enter  a  ray  convergence  routine.  The  horizontal 
distance  between  where  a  ray  Intersects  the  receiver  height  and 
the  receiver  range  is  given  by 

e  *  x  -  Range  (79) 

A  new  ray  is  traced  with  the  starting  angle 

8'  -  ©  -e  /Ox/38  )  (80) 

0  0  0 

where  8x/38q  is  given  by  equation  (72).  This  process  will,  under 
favorable  conditions,  reduce  the  value  of  £  ,  and  is  repeated  until 

e  is  smaller  than  a  specified  tolerance. 

B#  Some  examples 

The  present  models  may  be  used  to  analyze  a  multitude  of 
situations.  Only  a  few  can  be  discussed  here. 

First  to  be  considered  is  a  raised  maximum  phase  velocity. 

It  was  shown  in  section  II-C-7  that  this  would  cause  a  shadow 
zone.  The  question  discussed  here  is  how  intense  must  that 
maximum  be  to  show  a  noticeable  effect  and  also,  what  happens 
nearer  the  ground,  below  the  maximum,  since  rays  will  be  bent 
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downward.  In  figure  8,  there  is  an  iso-velocity  situation 
near  the  ground,  and  the  velocity  is  maximum  there.  In  the 
upper  portion  rays  are  bent  upward,  as  would  be  expected, 
toward  the  lower  velocity.  The  rays  contained  in  the 
iso-velocity  layer  are  straight  and  easily  penetrate  into 
the  upper  layer.  Figure  9  shows  a  slight  inversion  in  the 
lower  level.  The  same  rays  are  plotted  here  and  the  plot 
shows  that  the  rays  don't  penetrate  to  the  upper  layer 
quite  as  easily  as  before.  Rays  are  bent  downward  and 
trapped  by  the  inversion.  Figure  iO  shows  a  more  intense 
inversion.  The  rays  are  bent  as  before  but  to  the  right 
of  the  plot  are  more  concentrated  in  the  lower  part. 

Figure  11  shows  this  concentration  more  clearly.  More 
rays  have  been  added  between  the  rays  in  the  lower 
portion  of  figure  10.  It  is  noted  in  this  figure 
that  the  upper  section  is  much  more  concentrated  than 
the  lower,  indicating  a  much  higher  intensity  of  sound. 

This  point  may  be  considered  part  of  a  caustic.  It  is 
easily  seen  from  this  set  of  figures  that  the  more 
intense  an  inversion,  the  more  rays  may  be  trapped 
below.  This  would  indicate  that  the  sound  Intensity 
might  likely  be  much  higher  in  this  region. 

The  ray  tracing  program  may  be  used  with  a  variable 
terrain  as  seen  in  figure  12.  The  elgenray  routine  is 
not  yet  capable  of  this.  The  problem  is  that  the 
techniques  us  -o  to  speed  up  the  computation  time  take 
advantage  of  the  cyclic  nature  of  rays  which  exists 
only  if  there  are  similar  conditions  over  the  entire 
terrain.  Further  investigation  is  necessary  to  allow 
for  the  ability  to  handle  variable  topography  and 
maintain  optimal  use. 

Table  1  shows  the  output  of  the  elgenray  routine 
for  an  inversion  condition.  This  is  a  list  of  the  rays 
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TABLE  OF  EIGENRAYS 


TRAVEL 

TIME 

(SEC> 

START 

ANGLE 

(DEG) 

GROUND 
REFLECTIONS 
NO.  ANGLE 

ATTN 

LOSS 

(DB) 

SPREADING 

LOSS 

(DB) 

TOTAL 

LOSS 

(DB) 

O.Q 

5*117 

2 

9.120 

5.38 

68.19 

73.57 

0.0000 

-5.117 

2 

9.120 

5.38 

68.19 

73.57 

0.0058 

4.327 

3 

8.704 

5.39 

58.73 

64.11 

0.0058 

-4.327 

3 

8.704 

5.39 

58.73 

64.11 

0.0074 

2.807 

4 

8.061 

5.38 

58.55 

63.93 

0.0074 

-2.807 

4 

8.061 

5.38 

58.55 

63.93 

0.0083 

-1.996 

5 

7.817 

5.38 

57.72 

63.10 

0.0083 

1.996 

5 

7.817 

5.38 

57.72 

63.10 

0.0090 

1.406 

6 

7.688 

5.38 

56.50 

61.89 

0.0090 

-1.406 

6 

7.688 

5.38 

56.50 

61.89 

0.0097 

-0.916 

7 

7.614 

5.38 

54.42 

59.80 

0.0097 

0.916 

7 

7.614 

5.38 

54.42 

59.80 

0.0102 

0.391 

8 

7.569 

5.38 

48.52 

53.90 

0.0102 

-0.391 

8 

7.569 

5.38 

48.52 

53.90 

0.0825 

5.143 

l 

9.135 

5.38 

67.81 

73.20 

0.0940 

4.520 

2 

8.801 

5.39 

59.12 

64.51 

0.1228 

-0.585 

8 

7.582 

5.38 

50.26 

55.64 

0.1231 

2.978 

3 

8.122 

5.38 

58.99 

64.37 

0.1463 

2.171 

4 

7.863 

5.38 

58.35 

63.73 

0.1696 

1.592 

5 

7.724 

5.38 

57.46 

62.84 

0.1952 

1.135 

6 

7.644 

5.38 

56.15 

61.53 

0.2277 

0.728 

7 

7.594 

5.38 

53.86 

59.24 

0.2370 

-1.181 

7 

7.651 

5.38 

55.07 

60.45 

0.2835 

0.266 

8 

7.564 

5.38 

47.29 

52.67 

0.3557 

-1.800 

6 

7.770 

5.38 

56.91 

62.29 

0.5127 

-2.620 

5 

7.998 

5.38 

58.01 

63.40 

0.7611 

-4.117 

4 

8.602 

5.39 

58.25 

63.63 

1 . 2006 

-5.095 

3 

9.108 

5.39 

68.62 

74.00 

Table  1  -  List  of  elgenrays 
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that  Intersect  the  ease  receiver  point  specified  as 
nine-hundred  and  fifteen  aeters.  Figure  13  la  a  plot 
of  a  nuaber  of  these  rays  (froa  the  ray  tracing  prograa) 
and  shows  that  In  fact,  they  do  Intersect  at  the 
specified  receiver  location. 


' !  *  -w. 


Figure  13  -  Plot  of  elgenrays 
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IV.  Summary 

The  purpose  of  this  paper  has  been  to  discuss  ray-tracing 
techniques.  The  equations  have  been  derived  from  basic 
principles  in  a  straight  foward  manner.  Ray-tracing  may  be 
used  in  noise  control  applications  as  well  as  sound  ranging. 
Ray  solutions  are  good  approximations  to  wave  solutions 
under  the  condition  that  the  velocity  gradient  doesn't 
change  very  much  over  a  wavelength. 

An  analysis  of  the  ray  solution  has  been  performed. 
Caustics  are  formed  when  rays  are  either  bent  toward  each 
other  or  wavefronts  have  a  concave  profile.  Linear  theory 
predicts  that  there  is  infinite  energy  at  a  caustic.  This 
is  not  so  in  reality  due  to  non-linear  effects.  Caution 
must  be  taken  when  reviewing  output  from  a  ray  analysis. 
Although  the  theory  may  predict  infinite  energy  at  a 
caustic,  experiments  show  that  the  amount  of  energy  may 
be  very  large,  but  not  Infinite. 

Shadow  zones  occur  when  there  exists  an  effective 
maximum  sound  velocity  at  some  height.  Waveguides 
occur  when  there  is  a  raised  minimum  sound  velocity. 

In  anisotropic  media  rays  are  not  orthogonal  to 
wavefronts.  For  the  present  models  only  isotropic 
media  are  considered.  An  understanding  of  how  ray9 
travel  in  anisotropic  media  is  enlightening  to  real 
situations. 

Computer  programs  have  been  developed  to  demonstrate 
ray  techniques  and  are  contained  in  the  appendices  of  this 
paper.  These  programs  have  been  used  to  show  some  examples. 

The  programs  are  presently  being  utilized  in  much 
research  at  the  Moise  Control  Lab  of  The  Pennsylavnia 
State  University.  They  are  being  constantly  revised 
for  various  uses. 
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C  RAY  PATH  CALCULATION  -  MAIN  PROGRAM 

COMMON  /SX/DEP( 100) ,VEL( 100) ,GRAD<99) ,TEMP( 100) ,WV( 100) 
COMMON  /R/TT(99) , DB( 99 ) ,  ATN( 99) ,ANG0(99) ,ANGS(99) , ANGB ( 99 ) 
COMMON  /R/NS(99) ,NB(99) 

COMMON  /P/TLOSS( 99) 

INTEGER  TITLE 

DIMENSION  TITLE( 10) , BL( 10) 

5  READ(5 ,301 , END-6)  TITLE 

READ( 5 , 302 )  NDFT , NVFT , NTF , NWV , NRFT , I VEL , ITEMP , IWV , IRHM 
IF  ( ( I VEL+ITEMP ) . EQ . 0 )  STOP 

C  IF  SVP  DATA  IS  iO  BE  INTERNALLY  GENERATED,  REPLACE  'STOP'  BY 
C  APPROPRIATE  'GO  TO'  TO  GENERATING  ROUTINE. 

WRITE (6,305)  TITLE 
NP-0 

10  NP-NP+1 

READ ( 5 ,303)DEP(NP) ,TEMP(NP) ,VEL(NP) ,WV(NP),NOMO 
IF  (NOMO.EQ.O)  GO  TO  10 

15  CALL  SSP(NP, NDFT, NVFT, NTF, NWV, IVEL, ITEMP, IWV) 

READ (5 ,400)  NBL , POR, ( BL(I) , I- 1 , NBL ) 

READ (5,304)  TWIN 

20  READ( 5,304)  SD , TD , RANGE , ANGMAX .ANGMIN , FREQ , RHM 
C  IF  DEPTH  IN  FEET 

IF(NDFT.EQ. 1)  GO  TO  61 
C  CONVERT  TO  FEET 

SD-SD*3 . 28 
TD-TD*  3.28 

61  IF  (RANGE. EQ.O)  GO  TO  5 
C  IS  RANGE  IN  FEET 

IF(NRFT.EQ. 1 )  GO  TO  22 
C  CONVERT  METERS  TO  KILOYARDS 

RANGE-RANGE* 1 .09333/1000. 

GO  TO  23 

C  CONVERT  FEET  TO  KYARDS 

22  RANGE-RANGE/3000. 

23  IF( IRHM.EQ. 1  )  GO  TO  21 

C  DEFAULT  VALUE  OF  RELATIVE  HUMIDITY 

RHM-50 . 

21  IF ( ITEMP .EQ.O)  GO  TO  24 

C  FIND  AVERAGE  TEMPERATURE  IN  DEGREES  C 

AVTP-0.0 

IF( NTF .EQ.O)  GO  TO  26 
DO  27  I-l.NP 

27  AVTP-(TEMP(I)-32.)*5./9.+AVTP 
CO  TO  28 

26  DO  29  I- 1 , NP 
29  AVTP-TEMP ( I )+AVTP 

28  AVTP-AVTP/NP 
GO  TO  31 

C  DEFAULT  AVERAGE  TEMP  20  DEGREES  C 

24  AVTP-20 . 0 
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31  CALL  RAY (NP,SD,TD, RANGE , ANGMAX , ANGM1N , FREQ , NRAY , RHM , AVTP ) 
IF(NRAY.EQ.O)  GO  TO  71 
CALL  TTORD(NRAY) 

TO-TT(l) 

DO  25  1*1, NRAY 
25  TT( I)-TT( I)-TO 

SS-20.*ALOG10( 1E3*RANGE) 

C  IF  OUTPUT  IN  MRS  OR  BES 

IF(NVFT.EQ.  1 )  GO  TO  55 
C  CONVERT  TO  MRS 

SD-SD*. 304878 
TD-TD*. 304878 
RANGE-RANGE*. 9 146341 

WRITE(6 , 356 )  TITLE, SD.TD, RANGE, FREQ, ANGMAX .ANGMIN , SS , TO 
GO  TO  56 

55  WRITE( 6 , 353 )  TITLE , SD , TD , RANGE , FREQ , ANGMAX , ANGMIN , SS , TO 

56  DO  45  R-l.NBL 
DO  30  1*1, NRAY 
TLOSS ( I )-DB( I ) 

IF  (NB(I).EQ.O)  GO  TO  30 
XNB-NB ( I ) 

TLOSS(I)-DB(I)+XNB*BLOS(FREQ,POR,ANGB(I) )+XNB*BL(R)+ATN( I ) 
30  CONTINUE 

WRITE (6,450)  BL(R) 

WRITE(6 ,354) 

WRITE (6, 355)  (TT( I ) , ANGO (I) ,NB(I) ,ANGB(I) , 

1  ATN( I) , DB( I ) ,TLOSS(I) ,1*1 , NRAY ) 

CALL  IN TOUT (NRAY .TWIN ,XIOM) 

45  CONTINUE 

71  IF(NRAY.EQ.O)  WRITE(6,358) 

6  STOP 

400  FORMAT (II ,4X,11F5.1) 

450  FORMAT ( 1 5H  BOTTOM  LOSS  -  ,F5.1,//) 

301  FORMAT ( 1 0A4 ) 

302  FORMAT( 911) 

303  FORMAT( 4F 1 0.4, IX, II) 

304  FORMAT (7F10.3) 

305  FORMAT ( 1 H 1 , 1 0A4 / / ) 

353  FORMAT ( 1H1 , 10 A4// 

1  1H  , 1 2HS0URCE  DEPTH , F8 . 3 , 3H  FT/ 

2  1H  , 1 2HTARGET  DEPTH , F8 . 3 , 3 H  FT/ 

3  1H  , 5HRANGE , F8 . 3 , 5H  RYDS// 

4  1H  , 4HFREQ , F7 . 3 , 4H  RHZ / 

5  1H  , 9HMAX  ANGLE ,F6. I ,4H  DEG/ 

6  1H  , 9HMIN  ANGLE, F6. 1 ,4H  DEC// 

7  1H  , 1 2HSPH  SPP  L0SS,F7.2,3H  DB/ 

8  1H  ,  16H1 ST  ARRIVAL  TIME,F8.3,5H  SECS//) 

354  FORMAT (/// , 16X, 1 8 HT ABLE  OF  EICENRAYS// 

1  1H  , 20HTRAVEL  START  , 

2  35HGROUND  ATTN  SPREADING  TOTAL/ 


n  o 


i 
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3  1H  ,1811  TIME  ANGLE  , 

4  36HRE FLECTIONS  LOSS  LOSS  LOSS/ 

5  1H  , 1 9H  (SEC)  (DEG)  , 

6  35HNO.  ANGLE  (DB)  (DB)  (DB)//) 

355  FORMAT ( 1 H  , F6 . 4 , F9 . 3 , 2X , 14 , F8 . 3 , 2F8 . 2 , F  1 0 . 2 / ) 

356  FORMAT ( 1  HI , 10A4/ / 

1  1H  , 1 2HSOURCE  DEPTH , F8 . 3 , 3H  M  / 

2  1H  , 1 2HTARGET  DEPTH , F8 . 3 , 3 H  M  / 

3  1H  , 5HRANGE ,F8.3,5H  KM  // 

4  1H  , 4HFREQ , F7 . 3 , 4H  KHZ/ 

5  1H  , 9HMAX  ANGLE, F6. 1 ,4H  DEG/ 

6  1H  , 9HMIN  ANGLE, F6.1.4H  DEG// 

7  1H  , 1 2HSPH  SPP  LOSS , F7 . 2 , 3H  DB/ 

8  1H  , 1 6H 1  ST  ARRIVAL  TIME,F8.3,5H  SECS//) 

358  FORMAT ( 1  OX , 'NO  RAYS  FOUND') 

END 

SUBROUTINE  RAY ( NP , SD , TD , RANGE , ANGMAX , ANGMIN , FREQ , NRAY , RHM , AVTP ) 

C  PROGRAM  FINDS  EIGENRAYS  AND  CALCULATES  TRANSMISSION  LOSS 

C  NP  -  NUMBER  OF  POINTS  IN  SOUND  SPEED  PROFILE 

C  SD  -  SOURCE  DEPTH  (FT) 

C  TD  -  TARGET  DEPTH 

C  RANGE  -  SOURCE-  TARGET  HORIZONTAL  RANGE  (KYDS) 

C  ANGMAX  -  MAX  ANGLE  SEARCHED  (DEG) 

C  ANGMIN  -  MINUMUN  ANGLE  SEARCHED 

C  NRAY  -  NUMBER  OF  EIGENRAYS  FOUND 

C  RHM  -  RELATIVE  HUMIDITY  N  PER  CENT 
C  AVTP  -  AVERAGE  TEMPERATURE  N  DEGREES  CELSIUS 
C  AUX  PRINT-OUT:  SW7  ON  -  RAY  SEARCH  INFO 
C  SW8  ON  -  DF-BUG 

COMMON  /SX/D1( 100) ,V1(100) ,G1 (99) ,T1 1( 100) ,WV( 100) 

COMMON  /R/TT(99) ,DB( 99) , ATN (99) ,AN0(99) ,ANS( 99) , ANB( 99 ) 

COMMON  /R/LS(99) ,LB(99) 

DIMENSION  D( 102) ,V( 102) ,G(101) 

DIMENSIONDD(2) ,ND(2) 

DOUBLE  PRECISION  P ID , VKD , CVD , THOD , S ITHD , CSTHD , S ITH2 D , CSTH2D 

DOUBLE  PRECISION  XD , DXD ,XTD , RYARDD , SUMD  ,  DSUMD 

IPDB-2 

IPRINT-2 

PID-3. 14159265358979D0 
PI-SNGL(PID) 

C  MAX  NUMBER  OF  RAYS  (SIZE  OF  /R/  ARRAYS) 

NRAYMX-99 

CALCULATE  ATTN  COEFF  BY  AMERICAN  NATIONAL  STANDARD 
CHANGE  TO  DEGREES  KELVIN 
AVTP-AVTP+273. 15 
C  CHANGE  TO  HZ 

FTT-FREQMOOO. 

TO-293. 1  5 
TOl-273.  16 

PLR-10. 7 9 58 6* (  1  . -( TO  1 / AVTP ) )-5 . 0 2808*AL0G 1 0 ( AVTP /TO  1 )  +  l .  504  74*10.* 
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1  * ( -4 . )*(1.-10.**(-8.29692*((AVTP/T01)-1.)) )+0 .42873*10. **(-3. )*(-l 
l.+10.**(4. 7  6955  * (  1  .-(T01/AVTP) ) ) )-2. 2 19  5983 
HM-RHM*10.**PLR 

FRO-24.+4 ,41*10.**4*HM*(0.05+UM)/(0.391+HM) 

FRN-(T0/AVTP)**.5*(9.+350.*HM*EXP(-6. 1  42* ( ( AVTP/ TO ) ** ( - l ./3. )-l . )  ) 

1  ) 

C  ALPHA  IN  NEPERS/METER 

ALPHA-FTT**2*( 1 . 84* 1 0 . ** ( - 1 1 )*( AVTP/TO )**. 5+( AVTP/TO )**( -5 . /2.  )*( 1 
1 . 278*10.**(-2)*EXP(-2239. 1/AVTP)/ ( FRO+FTT* FTT/ FRO )+. 1068*EXP(-3352 
I . / AVTP ) / ( FRN+FTT*FTT/ FRN ) ) ) 

C  CONVERT  TO  DB/KYD 

ALPHA-ALPHA*868. 589*3.048037*3. 

C  FIT  SOURCE  AND  TARGET  INTO  SVP 
DO  5  J-1,NP 
D(J)-D1(J) 

V(J)-V1(J) 

IF(J.EQ.NP)  GO  TO  6 

5  C(J)«C1(J) 

6  LP-NP 
I-I 

IF  (SD-TD)  10,11,12 

10  DD(1)«SD 
DD(2)-TD 
J-l 

GO  TO  15  | 

11  DD(2)-SD  ' 

J-2 

GO  TO  15  l 

12  DD(1)-TD 
DD(2)-SD 
J-l 

15  IF  (DD(J)-D(I))  20,23,24 

20  LP-LP+1 
IP-LP-I 

DO  21  K-l.IP 
L-LP-K 
M-L+l 
D(M)-D(L) 

V(M)-V(L) 

IF(L.EQ.l)  GO  TO  26 
M-L-l 

21  G(L)-G(M) 

26  D(I)-DD(J) 

M-I-l 

V  (  I  )  -  V  (  M  )+G  (  M  )  *  (  D  (  I  )  -D  ( M  )  ) 

ND ( J ) -I 

22  IF  (J.GE.2)  GO  TO  35 
J-2 

GO  TO  15 

23  H D ( J ) - 1 
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GO  TO  22 

24  IF  (I.GE.LP)  GO  TO  30 
I-I+l 
GO  TO  15 
30  ND(J)-LP 

35  IF  (SD-TD)  40,41,42 

40  NSD-ND(l) 

NTD-ND(2) 

GO  TO  60 

41  NSD-ND(2) 

NTD-NSD 
GO  TO  60 

42  NTD-ND(l) 

NSD-ND( 2 ) 

C  INITIALIZE  RAY  TRACE 
60  ANGO-ANGMAX 

RYARD-1E3*RANGE 

RYARDD-DBLE(RYARD) 

ERRMX-l . 

STEP-0.05 

NSTEP-0 

IRAY-0 

J RAY-0 

NRAY-0 

IJ-0 

C  IPRINT-1  IF  SS7  ON:  IPRINT-2  IF  SS7  OFF 
IF  ( IPRINT.EQ.2)  GO  TO  65 

WRITE (6, 802)  SD,TD, RANGE , ANGMAX , ANGMIN , FREQ , ALPHA 
WRITE(6 ,801  ) 

IP-LP-1 

WRITE (6 ,800)  (I,D(I),V(I),G(I),I-1,IP) 

WRITE(6 ,800)  LP,D(LP) ,V(LP) 

WRITE(6 ,950) 

C  START  NEW  RAY 
65  K-NSD 

C  CHECK  IF  INITIAL  RAY  AT  HIGHEST  LIMIT 
IF  (K.GT.l)  GO  TO  70 

C  DOES  INITIAL  'LIMIT'  RAY  GO  DOWNWARD  ? 

IF  (ANGO.GT.O.)  GO  TO  205 

IF  ( (ANGO.EQ.O. ) .AND. (G( 1 ) .GE.O. ) )  GO  TO  205 
C  CHECK  IF  INITIAL  RAY  ON  GROUND 
70  IF  (K.LT.LP)  CO  TO  75 
C  DOES  INITIAL  GROUND  RAY  GO  UPWARD  ? 

IF  (ANCO.LT.O.)  GO  TO  210 

IF  ( (ANCO.EQ.O. ) .AND. (G(LP-1 ) .LE.O. )  )  GO  TO  210 
C  IS  INITIAL  ANGLE  ZERO  ? 

75  IF(ABS(ANGO) .GT. IE-3)  GO  TO  90 
C  IF  INITIAL  RAY  IS  SPLIT,  ARBITRARILY  MAKE  DOWNARD 
IF  ((C(K-l).CT.O.).AND.(G(K).LT.O.))  GO  TO  80 
C  IK  INITIAL  RAY  IS  DOWNARD  ,  DECREASE  ANGO  SLIGHTLY 
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IF  ( (G(K-1 ) .LE.O. ) .AND. (G(K) .LT.O. ) )  GO  TO  80 
C  IF  INITIAL  RAY  IS  UPWARD  ,  INCREASE  ANGO  SLIGHTLY 
IF  ( ( G (K-l ) . GT • 0. ) .AND .(G(K).GE.O*))  GO  TO  85 
C  MAKE  SPECIAL  CALCULATION  IF  RAY  IS  HORIZONTAL 
GO  TO  220 
80  AN GO ■ AN G 0-0 .01 
GO  TO  90 

85  ANC0-ANG0+0.01 
C  INITIALIZE  ANGO,  ETC 
90  THO-PI/ 180.*ANG0 
THOD-DBLE(THO) 

CSTHD-DCOS(THOD) 

CSTH-SNGL(CSTHD) 

CVD-DBLE(V(NSD)  )/CSTHD 
CV-SNGL(CVD) 

SITH-SIN(THO) 

SITHO-SITH 

X- 0.0 

XI- 0.0 
X2-0.0 
KV1-0 
KV2-0 

IBUG-  90;  IF(IPDB.EQ.l)  WRITE(6,888)  I  BUG , ANGO , S ITH , CSTH , CV 
C  CALCULATE  ONE  LAYER 

100  IBUG-100;  IF( IPDB.EQ . 1 )  WRITE(6,888)  I  BUG , V( K ) , SITH , SITH2 , X , X 1 , X2 
IF  (SITH. LT.O.)  GO  TO  110 
C  IF  RAY  GOES  UPWARD  BEYOND  LIMIT 

IF(K.LE.l)  GO  TO  205 
105  K-K-l 
DIR-1 . 

GRAD-G(K) 

GO  TO  120 

C  DOWNARD-GOING  RAY 

110  IF  (K.LT.LP)  GO  TO  115 
C  REFLECTION  OFF  GROUND 
SITH2--SITH2 

IF  (KV1.NE.0)  KV2-LP  ;  IF  (KV1.EQ.0)  KVl-LP 
GO  TO  140 
115  GRAD-G(K) 

K-K+l 

DIR--1 

C  DISTANCE  CALCULATION;  K  -  NEXT  LAYER 
C  ISO-VELOCITY  ? 

120  IF  (GRAD.EQ.O.)  GO  TO  125 
VKD-DBLE(V(K) ) 

IBUG-120;  IF( IPDB.EQ . 1 )  WRITE(6,888)  IBUC,V(K),CV 
C  VERTEX  IN  LAYER  K  ? 

IF  (VKD.GT.CVD)  CO  TO  130 
IF  (VKD.EQ.CVD)  GO  TO  205 
CSTH2*SNGL(VKD/CVD) 
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SITH2-DIR*SQRT( ( 1 . -CSTH2 ) * ( 1 • +CSTH2 ) ) 

DX-CV/GRAD* ( SITH2-SITH) 

GO  TO  135 

C  ISO-VELOCITY  CALCULATION 
125  ID-DIR 

LAST-K+ID 

SITH2-SITH 

CSTH2-CSTH 

DX-( D ( LAST ) -D ( K) )*CSTH2/SITH2 
GO  TO  135 

C  VERTEX  CALCULATION 
130  ID-DIR 
K-K-*-ID 

IF  (KV1.NE.0)  KV2-K  ;  IF  (KV1.EQ.0)  KVl-K 
SITH2--S ITH 
CSTH2-CSTH 
DX-2. *CV/GRAD*SITH2 
135  X-X+DX/3 
C  CHECK  RAY  POSITION 
C  RAY  AT  TARGET  DEPTH  ? 

IF  (K.NE.NTD)  GO  TO  140 
IF  (X1.GT.0.)  X2-X 
IF  (Xl.EQ.O.)  Xl-X 
C  RAY  RETURNED  TO  SOURCE  DEPTH  ? 

140  IBUG-140;  IF ( IPDB . EQ . 1 )  WRITE(6,888)  I  BUG , V ( K) , S ITH , S ITH2 , X , XI , X2 
IF  ( (K.EQ.NSD) .AND.(SITHO*SITH2.GT.O.  )  )  GO  TO  145 
IF( (X.GT. (1 .5*RYARD) ) . AND. (Xl.EQ.O.))  GO  TO  205 
SITH-SITH2 
CSTH-CSTH2 
GO  TO  100 
C  CYCLE  COMPLETED 
145  WL-X 

C  CHECK  1ST  INTERSECTION 

IF  (Xl.EQ.O.)  GO  TO  205 
NCYC-0 

ERRA-X1-RYARD 
150  ERRB-ERRA+WL 
C  MINIMUM  ERROR  NCYC  ? 

IF  (ABS(ERRB) .GE. ABS(ERRA) )  GO  TO  155 

ERRA-ERRB 

NCYC-NCYC+1 

IF  (NCYC. LT. 50)  GO  TO  150 
KIND-1 

IF  (IPRINT.EQ. 1 )  WRITE (6,902)  ANGO f KV 1 , KV2 , NCYC , KIND 
GO  TO  205 
C  1ST  RAY  ? 

155  IF  (IRAY.EQ.O)  GO  TO  160 
C  THIS  RAY  SAME  AS  LAST  ? 

IF  ((NCYC. EQ.ICYC). AND. (KV1.EQ.IV1). AND. (KV2.EQ.IV2))  GO  TO  170 
C  IF  NEW  RAY,  CALCULATE  INTENSITY  FOR  LAST  RAY 


►>  « - 
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CO  TO  280 
160  IRAY-IRAY+1 
ICYC-NCYC 
IVl-KVi 
IV2-KV2 

ERRIY-RYARD*1E60 
ERRIZ-ERRIY 
165  ANGI-ANGO 
ERRI-ERRA 
GO  TO  175 
170  ERRIX-ERRIY 
ERRIY-ERRIZ 
ERRIZ-ABS(ERRA) 

C  RANGE  ERROR  PASS  A  MAX  ? 

IF  ( (ERRIX.LT. ERRIY) .AND. (ERRIZ.LT. ERRIY)  )  GO  TO  280 
C  THIS  RAY  CLOSER  TO  TARGET  THAN  LAST  ? 

IF  (ABS(ERRA) .LT.ABS(ERRI) )  GO  TO  165 
C  CHECK  2ND  INTERSECTION 

175  IF  (X2.EQ.0.)  GO  TO  205 
NCYC-0 

ERRA-X2-RYARD 
180  ERRB-ERRA+WL 

IF  (ABS(ERRB) .GE.ABS(ERRA) )  GO  TO  185 

ERRA-ERRB 

NCYC-NCYC+1 

IF  (NCYC.LT.50)  GO  TO  180 
KIND-2 

IF  (IPRINT.EQ.l)  WRITE (5,902)  ANGO , KV 1 , KV2 , NCYC , KIND 
GO  TO  205 

185  IF( JRAY.EQ.O)  GO  TO  190 

IF  ((NCYC. EQ.JCYC). AND. (KV1.EQ.JV1). AND. (KV2.EQ.JV2))  GO  TO  200 
GO  TO  285 
190  JRAY-JRAY+1 
JCYC-NCYC 
JV1-KV1 
JV2-KV2 

ERRJY-RYARD* 1 E60 
ERRJZ-ERRJY 
195  ANGJ-ANGO 
ERRJ-ERRA 
GO  TO  205 
200  ERRJX-ERRJY 
ERRJY-ERRJZ 
ERRJZ-ABS ( ERRA) 

IF  ( (ERRJX.LT. ERRJY) .AND. (ERRJZ.LT. ERRJY) )  GO  TO  285 
IF  (ABS(ERRA) .LT.ABS(ERRJ) )  GO  TO  195 
C  DECREMENT  ANGO 
205  NSTEP-NSTEP+1 

STEPN-FLOAT(NSTEP) 

ANGO-ANCMAX-STEPN*STEP 
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C  DECREMENTED  THRU  THE  RANGE  ? 

IF  (ANGO.GE.ANGMIN)  GO  TO  65 
C  CONVERGE  LAST  I  AND  J  RAYS 
210  IJ-1 

C  IF  NO  RAYS  FOUND 

IF( (IRAY+JRAY+NRAY) .EQ.O)  RETURN 
GO  TO  280 

C  HORIZONTAL  RAY  CALCULATION 
220  IF  (NSD.NE.NTD)  GO  TO  205 
ERR»0. 

S-RANGE 

TIM-3. *RYARD/V(K) 

SPL-20.*ALOG10(RYARD) 

ATTN-ALPHA*S 

NS-0 

NB-0 

LCYC-0 

KIND-0 

LV1-K 

LV2-K 

WRITE (6, 951)  ANGO, ERR, NS, NB.S, TIM, SPL, ATTN, LV 1 .LV2.LCYC, KIND 
GO  TO  205 

C  ZERO  IN  ON  TARGET  AND  CALCULATE  INTENSITY  LOSS 
280  KIND-1 

ANGL-ANG I 
LCYC-ICYC 
ERRL-ERRI 
LV1-IV1 
LV2-IV2 
GO  TO  290 
285  KIND-2 

ANCL-ANGJ 

LCYC-JCYC 

ERRL-ERRJ 

LVl-JVl 

LV2-JV2 

IF  (IJ.EQ.l)  IJ-2 
290  THOD-PID*DBLF,(ANGL)/  180D0 
THO-SNGL(THOD) 

ERRP-2 . *ERRL 
IVTX-0 
295  NS-0 
NB-0 
INT-O 
MV  1-0 
MV2-0 
ANGB-0.0 
ANGS-0.0 
X-0.0 
XT-0.0 
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S-0.0 

TIM-0.0 

SUll-O.  0 

X0-000 

XTD-ODO 

SUMD-ODO 

SITHD-DSIN(THOD) 

CSTHD-DCOS(THOD) 

CVD-DBLE(V(NSD) )/CSTHD 
CV-SNGL(CVD) 

SITHO-SNGL(SITHD) 

CSTUO-SNGL(CSTHD) 

TNTHO-SITHO/CSTHO 

ANGLO- 180. /PI*THO 

TH-THO 

SITH-SITHO 

CSTH-CSTHO 

K-NSD 

300  IF  (SITH.LT.O.)  GO  TO  310 
IF  (K.LE.l)  GO  TO  205 
305  K-K-l 
DIR-1 . 

GRAD-G(K) 

GO  TO  320 

310  IF  (K.LT.LP)  GO  TO  315 
NB-1 

IF  (MV1.NE.0)  MV2-LP  ;  IF  (MV1.EQ.0)  MV1-LP 
S2-SITH*SITH 

ANGB- 180. /PI* ATAN (SQRT(S2/ ( 1 . -S2 )  )  ) 

SITH2D— SITH2D 
SITH2  — SITH2 
TH2— TH2 
GO  TO  355 


315  CRAD-C(K) 
K-K+l 
DIR— 1 . 


320  IF 

(GRAD.EQ.O.  ) 

GO 

TO 

335 

VKD 

-DBLE(V(K)  ) 

IF 

(VKD.GT.CVD) 

GO 

TO 

340 

IF 

(VKD.EQ.CVD) 

GO 

TO 

379 

CSTH2D-VKD/CVD 

SITH2D-DBLE(DIR)/CVD*DSQRT( ( CVD-VKD ) * ( C VD+VKD ) ) 
CSTH2-SNGL(CSTH2D) 

SITH2-SNGL(SITH2D) 

TH2-ATAN (SITH2/CSTH2) 

325  DXD-CVD* ( SITH2D-SITHD ) / DBL£( GRAD ) 

DX-SNGL ( DXD ) 

DS»CV/GRAD*(TH2-TH) 

ARG-SNGL( ( 1D0+SITH2D)/ ( 1D0-SITH2D)*( 1D0-SITHD)/ ( 1D0+SITHD) ) 
DTIM-O. 5/GRAD*AL0G( ARC) 
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330  DSUMD-DXD/SITH2D/SITHD 
XD-XD+DXD/3D0 
X-SNGL(XD) 

S-S+DS/3. 

TIM-TIM+DTIM 

SUMD-SUMD+DSUMD/3D0 

SUM-SNGL(SUMD) 

GO  TO  345 

335  ID-DIR 

LAST-K+ID 

TH2-TH 

SITH2-SITH 

CSTH2-CSTH 

SITH2D-SITHD 

CSTH2D-CSTHD 

H-D(LAST)-DOC) 

DXD-DBLE(H)*CSTH2D/SITH2D 

DX-SNGL(DXD) 

DS-SQRT(DX*DX+H*H) 

DTIM-DS/V(K) 

GO  TO  330 

340  ID-DIR 
K-K+ID 
TH2--TH 
SITH2--SITH 
CSTH2-CSTH 
SITH2D  — SITHD 
CSTH2D-CSTHD 

IF  (MVl.NE.O)  MV2-K  ;  IF  (MV1.EQ.0)  MV1-K 
GO  TO  325 

345  IF  (K.NE.NTD)  GO  TO  355 
INT-INT+1 

IF  ( INT.NE.KIND)  GO  TO  355 

XTD-XD 

XT-SNGL ( XTD ) 

ST-S 

TIMT-TIM 

SUMT-SUM 

NST-NS 

NBT-NB 

355  IF  ( (K.EQ.NSD) .AND* (TH0*TH2 . GT . 0 . ) )  GO  TO  360 

IF  (((X.GT.(l. 5*RYARD ) ) .AND.(INT.EQ.O) ) .OR. ( INT.GT. 2) )  GO  TO  375 

TH-TH2 

SITH-SITH2 

CSTH-CSTH2 

SITHD-SITH2D 

CSTHD-CSTH2D 

GO  TO  300 

360  CYCL-FLOAT(LCYC) 

XD-XTD+XD*DBLE(CYCL) 
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X-SNGL(XD) 

ERR*SNGL ( XD-RYARDD ) 

IF  ((MV1.NE.LV1).0R.(MV2.NE.LV2))  GO  TO  378 

IF  (ABS(ERR) .GE.ABS(ERRP) )  GO  TO  370 

SUM«SUMT+CYCL*SUM 

IF  (ABS(ERR) .LE.ERRMX)  GO  TO  365 

DXDTH— SUM*TNTHO 

DTHO-ERR/DXDTH 

DANGO* 1 80 . / PI *DTHO 

IF  ( ABS ( DAN GO ).GT.(10.*STEP))  GO  TO  377 
THOD-THOD-DBLE(DTHO) 

THO-SNGL(THOD) 

ERRP-ERR 
GO  TO  295 
365  S-ST+CYCL*S 

TIM“T IM+CYCL*TIM 

SPL«10.*ALOG10( ABS ( X*S ITH2*TNTH0*SUM/ CSTHO ) ) 

NS“NST+LCYC*NS 

NB-NBT+LCYC*NB 

S-IE-3*S 

ATTN-ALPHA*S 

IF  (IPRINT.EQ. 1)  WRITE(6,951) 

1  ANGLO, ERR, NS, NB.S, TIM, SPL, ATTN, LV 1 ,LV2 ,LCYC, KIND 
NRAY-NRAY+1 
TT(NRAY)-TIM 
DB(NRAY)-SPL 
ATN(NRAY)-ATTN 
ANO(NRAY) -ANGLO 
ANS ( NR AY )"ANGS 
ANB(NRAY)-ANGB 
LS ( NRAY )“NS 
LB(NRAY)-NB 

IF  (NRAY.LT.NRAYMX)  GO  TO  380 
WRITE(6 ,805) 

RETURN 

370  IF  ( IPRINT.EQ. 1 )  WRITE(6,952)  ANGL , ANGLO , LV I , LV2 , LCYC , KIND 
GO  TO  380 

375  IF  ( IPRINT. EQ. 1 )  WRITE(6,953)  ANGL , ANGLO , LV I , LV2 , LCYC , KIND 
GO  TO  380 

377  IF  (  IPRINT. EQ.l)  WRITE(6,954)  ANGL , DANGO , LV 1 , LV2 , LCYC , KIND 
GO  TO  380 

378  IF  (IVTX.GE.3)  GO  TO  379 
IVTX-IVTX+I 
DTH0-DTH0/2. 

THOD»THOD+DBLE(DTHO) 

THO-SNGL(THOD) 

GO  TO  295 

379  IF  ( IPRINT.EQ. 1 )  WRITE(6,955)  ANGL , ANGLO , LV 1 , LV2 , LCYC  ,  KIND 

380  IF  ((IJ. EQ.l). AND. (JRAY.GT.O))  GO  TO  285 

IF  ( ( IJ . EQ • 2 ) .OR. ((IJ.EQ.l) .AND • ( JRAY . EQ . 0 ) ) )  RETURN 


IF  (KIND.EQ.l)  GO  TO  160 
GO  TO  190 

800  FORMAT ( 1 H  , 12 , 2F 1 0 . 2 , F 1 2 . 3  ) 

801  FORMAT ( 1 H  ,32HSVP  WITH  SOURCE  AND  TARGET  ADDED//) 

802  FORMAT( 1H1 , 12HSOURCE  DEPTH , F8 . 1 , 3H  FT/ 

1  IH  , 1 2HTARGET  DEPTH , F8 . 1 , 3H  FT/ 

2  IH  , 5HRANGE , F8 . 3 , 5H  KYDS// 

3  IH  , 9HMAX  ANGLE, F6. I ,4H  DEG/ 

4  IH  , 9HMIN  ANGLE, F6. 1 ,4H  DEG// 

5  IH  ,9HFREQUENCY,F7.3,4H  KHZ/ 

6  IH  , 1 OHATTN  COEFK , 1  PE  1 0 . 2 , 7H  DB/KYD////) 

805  FORMAT ( 1H0.48H***  FOUND  TOO  MANY  RAYS  -  DECREASE  ANGMAX , ANGMIN 
902  FORMAT ( 1 H  ,F7.3,21H  CYCLE  LIMIT  EXCEEDED , 3 2X , 1 5 , 1 7 , 1 7 , 16/ / ) 

950  FORMAT ( 1H1 , 17X, 37 HT ABLE  OF  SOUND  RAY  PATH  INTERSECTIONS// 

1  IH  ,3711  INITIAL  RANGE  NUMBER  OF  RAY  , 

2  44HTRAVEL  SPREADING  ATTN  1ST  2ND  NUMBER/ 

3  IH  , 38HANGLE  ERROR  REFLECTIONS  LENGTH  , 

4  47HTIME  LOSS  LOSS  VERTEX  VERTEX  OF  RAY/ 

5  IH  , 3 7 H  (DEG)  (YDS)  SURFACE  BOTTOM  (KYDS)  , 

6  49H( SECS )  (DB)  (DB)  LAYER  LAYER  CYCLES  TYPE 


) 


//) 


951  FORMAT ( 1 H  ,F7.3,F6.1,I6,I?,F9.2,F8.3,F8.2,F9.3,I5,I7,I7,I6//) 

952  FORMAT ( 1 H  ,F7.3,16H  RAY  DIVERGED  AT,F9.3,4H  DEG , 24X, 1 5 , 1 7 , 17 , 16/ / ) 

953  FORMAT ( 1 H  ,F7.3,12H  RAY  LOST  AT,F9.3,4H  DEG , 28X , 15 , 17 , 17 , 16 / / ) 

954  FORMAT'OH  ,  F7 . 3 , 1  5H  ATTEMPTED  JUMP,F9.3,4H  DEG,  25X,  15 , 17 , 17 , 16// ) 

955  FORMAT ( 1 H  ,F7.3,15H  DIFF  VERTEX  AT,F9.3,4H  DEG , 2 5X , 15 , 1 7 , 1 7 , 16 / / ) 
888  FORMAT (1H0,I8/10(1PE13.5)) 

END 

FUNCTION  BLOS(F,P, THETA) 

C  CALCULATES  BOTTOM  LOSS  FROM  NUWC  TECH  NOTE  10  (DEC  67). 

C  F  -  FREQ  (KHZ),  P  -  POROSITY,  THRTA  -  BOTTOM  CRAZING  ANGLE 

DIMENSION  ABTLOS (14) 

DATA  ABTLOS ( 1 ) , ABTLOS ( 2 ) .ABTLOS ( 3 ), ABTLOS ( 4 ), ABTLOS ( 5 ) , 

1  ABTLOS(6) ,ABTLOS(7) ,ABTL0S(8) ,ABTLOS(9) , ABTLOS ( 10) , ABTLOS ( 11), 

2  ABTLOS (12) , ABTLOS (13) , ABTLOS (  14)  /. 16, .67,1. ,1.18, 1.31, 1.43, 1.52, 

3  1.61,1.7,1.76,1.82,1.88,1.94,2./ 

BLOS-O.O 

IF(P.LT.O.Ol)  RETURN 
IF(F.GT.O.l)  GO  TO  15 
FUNU-0. 16 
GO  TO  40 

IF(F.LT.6.5)  GO  TO  20 
FUNU-2.0 
CO  TO  40 
DO  30  1-1,7 
XI-I 

IF(XI*0.5.GT.F)  GO  TO  35 
30  CONTINUE 

35  IF(I.EQ.l)  GO  TO  45 

FUNU-ABTLOS( I )+( ABTLOS ( 1+1 ) -ABTLOS ( I))*(F-(XI-1 . )*0.5)/0.5 
GO  TO  40 


15 


20 
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4  5  FUNU- ABTLOS ( 1 )  +  ( ABTLOS ( 2 ) - ABTLOS ( 1 ) ) * ( F-0 . 1 ) /O  .  4 
40  ARG- 1 . 5/ P*  ALOG ( P* THETA/ 13.74) 

ARG-EXP(ARG) 

BL0S-(3 .7+17.5*(P-.27))*FUNU*(TANH(ARG)+( 1 . 0-P/0 . 27 ) / 1 2 . 5* 
1  (THETA/90. 0)**2) 

RETURN 

END 

SUBROUTINE  INTOUT ( NRAY .TWIN , XIOM ) 

C  SUMS  INTENSITY  IN  MOVING  WINDOW  TWIN  SECONDS  LONG. 

COMMON  /R/TT(99) ,DB(99) ,ATN(99) ,ANGO(99) ,ANGS(99) , ANGB( 99 ) 
COMMON  /R/NS (99 ) , NB( 99 ) 

COMMON  /P/TLOSS(99 ) 

DIMENSION  XINT( 99 ) 

XLN10-0. 23025851 
XIOM— 400. 

WRITE ( 6 ,400  )  TWIN 
DO  10  I-l.NRAY 

10  XINT(I)-EXP(-XLN10*TL0SS(I)) 

SUM-0. 

Kl-1 
K2- 1 

15  IF(TT(K1 )-TT(K2)+TWIN)  30,30,20 
20  T2-TT ( K2 ) 

SUM-SUM+XINT ( K2 ) 

L2-K2 
K2-K2+1 
GO  TO  40 

30  T2-TT(K1 )+TWIN 
SUM-SUM-XINT(K1 ) 

Ll-Kl 

Kl-Kl+1 

IF( (TT(L1)-TT(X2)+TWIN) .EQ.O.  )  GO  TO  20 
40  RCV- 1 0 . * ALOG 10 ( ABS ( SUM+1 E-30 ) ) 

WRITE (6,401)  T2 , K 1 , L2 , RCV 
XIOM- AMAX1( XIOM, RCV) 

IF(K2 .LE.NRAY)  GO  TO  15 
WRITE (6,450)  XIOM 

450  FORMAT ( 1 H  , 24HMAX  INTEGRATOR  OUTPUT  -  ,  F  1 0 . 2  ,///// ) 

RETURN 

400  FORMAT( III , 2X, 1 7 H INTEGRATOR  OUTPUT/// 

1  1H  , 1 1HTIME  WINDOW, F6. 3, 4H  SEC// 

2  1H  , 23H  TIME  1ST  2ND  OUTPUT/ 

3  1H  , 22H  (SEC)  RAY  RAY  (DB )//) 

401  FORMAT ( 1 H  , F6 . 4 , 14 , 14 , F9 . 2  ) 

END 

SUBROUTINE  TTORD(NRAY) 

C  ORDERS  EIGENRAYS  BY  TRAVEL  TIME 

COMMON  /R/TT(99),DB(99),ATN(99) , ANGO ( 99 ) , ANGS ( 99 ) , ANGB ( 99 ) 
COMMON  /R/NS(99) ,NB(99) 

IE-NRAY-1 


o  u  o  o  u  o 
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DO  23  I-l.IE 
JS-I+1 

DO  25  J-JS.NRAY 
IF  (TT(J)  .GE.TT(I))  GO  TO  25 
TEMP  l-TT(l) 

TEMP  2-DB(I) 

TEMP  3-ATN(I) 

TEMP  4-ANG0(I) 

TEMP  5-ANGS(I) 

TEMP  6-ANGBd) 

NTEMP  l-NS(I) 

NTEMP  2-NB(I) 

TT(I)-TT(J) 

DB(I)-DB(J) 

ATM(I)-ATN(J) 

ANG0(I )-ANG0( J) 

ANGSd)-ANGS(J) 

ANGB(I)-ANGB( J) 

NS(I)-NS(J) 

NB(I)-M3(J) 

TT ( J ) “TEMP  1 
DB( J)-TEMP2 
ATN(J)-TEMP3 
ANGO( JJ-TEMP4 
ANGS ( J ) “TEMP 5 
ANGB( J)-TEMP6 
NS ( J ) "NTEMP 1 
NB( JJ-NTEMP2 
25  CONTINUE 
RETURN 
END 

SUBROUTINE  SS P ( NP , NDFT , NVFT , NTF , NWV , IVEL , ITEMP , IWV ) 
CALCULATE  SOUND  SPEED  PROFILE  FROM  BERANAK 
D-DEPTH 

G»SOUND  SPEED  GRADIENT 
V-SOUND  SPEED 
T-TEMP 

WV-WIND  VELOCITY 

COMMON / SX/ D(IOO ) ,V(100) ,G( 99) , T( 100 ) , WV (100 ) 

C  SOUND  SPEED  GIVEN  ? 

IF  (IVEL.EQ.l)  GO  TO  50 
C  TEMP  IN  DEG  F  ? 

IF  (NTF.NE.l)  GO  TO  10 
C  CONVERT  TEMP  TO  DEG  C 
DO  5  I-l.NP 

5  T(I)-5 ./9 .*(T(I)-32.  ) 

C  DEPTH  IN  FT  ? 

10  IF  (NDFT.NE.l)  GO  TO  20 
C  CONVERT  DEPTH  TO  METERS 
DO  15  1*1 , NP 
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15  D<I)-D(I)*0.3048 
C  WIND  VELOCITY  GIVEN? 

20  IF(IWV.EQ.l)  GO  TO  27 
C  ASSUME  NO  WIND  VELOCITY 
DO  25  I-l.NP 
25  WV(I)-0.0 

C  WIND  VELOCITY  IN  METERS/SEC? 

27  IF(NWV.EQ.O)  GO  TO  30 
C  CONVERT  TO  METERS/SEC 

DO  2  6  I* i  ,NP 

26  WV(I)-WV(I)*.  304878 

C  CALCULATE  SOUND  SPEED 
30  DO  35  1*1, NP 

35  V( I)-331 .4*(SQRT(1 .0+0. 00366 *T(I) ) )+WV(I) 
C  CALCULATE  SOUND  SPEED  GRADIENT 
DO  40  1-2, NP 
J-l-  1 


40  G(J)-(V(I)-V(J)  )/ (D(I)-O(J)  ) 

C  IF  VELOCITY  INPUT  IN  MRS  OUTPUT  IN  MRS 
IF(NVFT.EQ.O)  GO  TO  45 

46  DO  47  I-l.NP 
V(I)-V(I)*3. 28084 
D(I)«D(I)*3. 28084 
WV(I)-WV{ I)*3. 28084 

47  T(I)-9./5.*T(I)+32. 

IF(NVFT.EQ.O)  RETURN 

C  PRINT  PROFILE  INCLUDING  TEMP  AND  WV 
WRITE (6 ,801) 

IF-NP-1 
WRITE(6 ,802) 

WRITS(6 ,802) 

RETURN 

C  MRS  OUTPUT 

45  WRITE(6 ,805) 

IF-NP-1 
WRITE(6 ,802) 

WRITE(6 ,802) 

GO  TO  46 

C  SOUND  SPEED  IN  FPS 
50  IF  (NVFT.EQ.l) 


(I,D(I) ,V(I) ,T(I) ,  WV( I ) ,G(I) ,1-1 ,IF) 
NP,D(NP) ,V(NP) ,T(NP) ,WV(NP) 


(I,D(I),V(I),T(I),WV(I),G(I),I-l,IF) 
NP,D(NP) ,V(NP) ,T(NP) ,WV(NP ) 


GO  TO  60 


C  OUTPUT  IN  MRS  UNITS 

IF(NDFT)  53,70,53 
53  DO  51  I-l.NP 

51  D(I)-D(I)*.  304878 
NDFT-0 

GO  TO  70 

C  CONVERT  SPEED  TO  FPS 

52  DO  55  I-l.NP 

55  V(l)-V(I)*3. 28084 
C  DEPTH  IN  FT  ? 
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60  IF  (NDFT.EQ.l)  GO  TO  70 
C  CONVERT  TO  FEET 
DO  65  1*1 ,  NP 
65  D(I)-D(I)*3. 28084 
IF (NVFT.EQ.O) RETURN 
C  CALCULATE  SOUND  SPEED  GRADIENT 
70  DO  75  1-2, NP 
J-I-l 

75  G(J)*(V(I)-V(J))/ (D(I)-D(J)) 

C  PRINT  PROFILE 
C  MRS  OUTPUT  OR  3ES 

IF(NVFT.EQ. 1 )  GO  TO  80 
WRITE(6 ,806) 

IE-NP-i 

WRITE (6 ,804)  (I ,D(I) ,V(I) ,G(I)  ,1-1  ,IE) 

WRITE(6 ,804  )  NP,D(NP),V(NP) 

GO  TO  52 

80  WRITE (6 , 803  ) 

IE-NP-1 

WRITE <6 ,804)  (I,D(I),V(I),G(I),I-l,IE) 

WRI TE ( 6 , 804 )  NP ,D(NP) , V(NP) 

RETURN 

801  F0RMAT(20X, 19HSOUND  SPEED  PROFILE// 

1  8X, 4 1 HDEPTH  SPEED  GRADIENT  TEMP  WV/ 

2  9X , 44H( FT )  (FT/SEC)  (FPS/FT)  (DEG  F)  (FT/SEC)//) 

802  FORMAT ( 1 H  , 1 2 , 1  PE  1 0 . 3 , OP  FI  0 . 3 , 1 1 X , 2F8 , 3/ 2 1 X , 1  PE  1 1 . 3 ) 

803  FORMAT( 10X, 19HSOUND  SPEED  PROFILE// 

1  8X , 2  6HDEPTH  SPEED  GRADIENT/ 

2  9X, 25H( FT)  (FT/SEC)  (FPS/FT)//) 

804  FORMAT ( 1  H  , 12 , OPE  1 0 . 3 , F 10 . 3 / 2 1 X , 1  PE  1 1 . 3 ) 

803  FORMAT(20X, 19HS0UND  SPEED  PROFILE// 

1  3X , 4 1 HDEPTH  SPEED  GRADIENT  TEMP  WV/ 

2  9X,42H(M)  (M/SEC)  (M/S/M)  (DEG  C)  (M/SEC)//) 
806  F0RMATC10X, 19HS0UND  SPEED  PROFILE// 

1  8X , 26HDEPTH  SPEED  GRADIENT/ 

2  9X ,  2  5H (M)  (M/SEC)  (M/S/M)  //) 

END 


n  on  on 


Appendix  B 


C  SRP-PLOT  2  PROGRAM 

C  (4/6/73) 

C  COMPUTES  SOUND  RAY  PATHS  AND  WRITES  TAPE  OR  PUNCHES  CARDS  FOR 
C  PLOTTING  PATHS  USING  CALCOMP. 

C  WILL  TAKE  UP  TO  150  BOTTOM  DEPTH  COORDINATES  AND  5  SVPS  WITH  A 

C  MAXIMUM  OF  200  POINTS  IN  A  SVP.  EACH  SVP  MUST  REACH  THE  SAME 

C  MAXIMUM  DEPTH  AS  THE  ADJACENT  GIVEN  SVPS. 

DOUBLE . PRECISION  ALPH , ANGLE ( 3000 ) , ANGX , AQ , B , BDEP ( 1 50 ) , BETA, BQ, 

1  BRAN ( 150) , CQ ,CSTH , CSTHX, CSTH2 , D(6 , 1 00 ) , DD , DDEPTH , DELR , DELX , DELY , 

2  DEP ( 3000 ),DEPMIN,DISC,DM,DR,DRFT,G(6,100) , GAM , GC , GI , GI 2 , GRAD , P , 

3  P2,P3,P4,R(3000) ,RA,RC,RC2 , RMAX , RMAXL , RS VP ( 6 ) ,SAL(6, 100) ,SC,SD, 

4  SITH.SITHX, SITH2 .SLOPE ( 1 50 ) , SR , T ( 6 , 1 00 ) , TATH , TANTH , TANTHX , TC , TC2 

5  TC3 ,TC4 ,TD ,TEMP,TEM2 ,TEM3,TEM4 ,TF(6 , 100) , THONE , TR , V ( 6 , 100) , VC, VP 

6  VS,VSTP,VT,X, XBRAN , XP , Y , WV ( 6 , 100) 

INTEGER* 2  MA( 1000) ,MB( 1000) ,MC( 1000 ) 

INTEGER  IGSVP(6,250) ,NPSVP(6) , GRAPH/ 0/ 

REAL  DQ ( 6 , 1000) ,QD(6) , PI/3. 141593/ 

LOGICAL* 1  TITLE(40) , ALPD (20),TSVP(5) 

FUNC(A,B,C,D,E)  -  A  +  (B  -  C)*(E  -  A)/(D  -  C) 

DATA  TSVP/  '  S' , 'V'  , , '  ' / 

CALL  NOPRQ 

CALL  INITQ(MA,MB ,MC , DQ , 1000 ) 

10  READ (5,701)  TITLE 
HTMIN-0 .00 
HTMAX* 100.0 


READ ( 5 , / 0  2 )  NRSVPS ,NRBOT, METER, NAUT , DELR, NOUT, NEW SVP , NOP R, NO AD, 
1  NRAN , RANINC , RANL , NDEP , DEP INC , DEPL , NS  V , S V  INC , S VL , SVM IN, RANGE 
IF ( NOUT . LT  .  7  1  )  GO  TO  14 
CALL  STSWQ (564,71  ) 

14  IF ( NOAD  .LT  .  1  )  GO  TO  16 
READ ( 5,701 )  ALPD 
16  IF(N EWSVP.LT. 1  )  GO  TO  133 
DO  128  NR  -  1, NRSVPS 

READ (5, 7 03)  NODEP , NODFT , NOTEMP , NOTF , NO VEL , NO VFT , NO SAL 
NS VP  -  1 

20  READ (5,704)  D ( NR , NS VP ) , TF ( NR , NS VP ) , WV ( NR , NS VP ) , S AL ( NR , NS V P ) , 

1  IGSVP(NR.NSVP) ,NOMO 
D(NR,NSVP ) “HTMAX- ( D(NR,NSVP) -HTM IN ) 

IF(NOMO.GT.O)  GO  TO  24 
NSVP  -  NSVP  +  1 
GO  TO  20 

24  DO  97  1*1 , NS  VP 

97  V ( NR , I ) ■ 3  3 1 . 4* ( DSQRT (  1  .0+0 .0036  6*TF( NR , I ) ) )+WV ( NR , I ) 

98  J  -  NSVP  -  1 

CALCULATE  VELOCITY  GRADIENTS 


DO  103  I  -  l.J 


o  o 
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G ( NR , I )  -  ( V (NR, 1+1 )  -  V(NR , I ) ) / ( D ( NR , 1+1 )  -  D(NR,I)) 

IF(DABSCG(NR, I ))  .GE. 0.0002)  GO  TO  103 

G  (  NR ,  I  )  -  0.0 
103  CONTINUE 

G(NR,NSVP)  -  G(NR.NSVP-l) 

IB  -  1 
NP AGE  -  l 

107  WRITE (6,705)  N PAGE, TITLE, NR, NSVP 
IF ( NOVEL . LT  .  1  )  GO  TO  110 
WRITE( 6,706) 

110  WRITE( 6,707)  NR 
IL  -  IB  +  21 
IF(IL.LE.J)  GO  TO  114 
IL  »  J 

114  IF(NOTEMP.GT.O)  GO  TO  117 

WRITE (6 , 708)  (D(NR, I) ,TF(NR, I ) , SAL (NR, I ) , V(NR, I) ,G(NR, I) , 

1  I  «  IB , IL) 

GO  TO  118 

117  WRITE (6,709)  ( D ( NR , I ) , V ( NR , I ) , G ( NR , I ) ,  I  *  IB,IL) 

'.18  IB  -  IL  +  1 

IF(IB.GT.J)  GO  TO  123 
WRITE(6 ,710) 

NPAGE  -  NP AGE  +  1 
GO  TO  107 

123  IF(NOTEMP.GT.O)  GO  TO  126 

WRITE (6, 708)  D( NR, NS VP) , TF (NR, NSVP), SAL (NR, NSVP), V(NR, NSVP) 

GO  TO  127 

126  WRIT E(6, 709)  D ( NR , NS VP ) , V ( NR , NS VP ) 

127  NPSVP(NR)  -  NSVP 

128  CONTINUE  j 

LABEL  PLOTS  ******************************** 

DO  132  I  -  1 , NS VP  | 

D(NRSVPS  +  1  ,  I)  -  D( NRSVPS  ,  I  ) 

G ( NRS VPS+l , I )  -  G( NRSVPS,  I) 

132  V ( NRSVPS+1 , I  )  -  V ( NRSVPS , I  )  , 

133  DISH  -  0.5MRANL  +  FLOAT  ( NRSVPS  )*(  SVL  +  1.0)) 

CALL  XFSTQ(DISH, 1.6, 0.24,0. 26, 0.0,0.0,QD) 

CALL  LABLQ(T1TLE,-40,OD,GRAPH,40) 

DISSR  -  RANL/2.0 

IF ( NOAD . LT . 1 )  GO  TO  143 

CALL  XFSTQ(  DISSR  -  1  .  7 , 0 . 9 5 , 0  .  1  8 , 0 . 2 0 , 0 . 0 , 0 . 0 , QD ) 

CALL  LABLQ( 'SOUND  RAY  PATHS  -  INTENSITY  CON  TO U RS '  , - 3 8 , Q D , G RAP H , 3 8 )  j 
CALL  XFSTQ(DISSR  +  3 . 7 , 0 . 9 5 , 0 . 1 5 , 0 . 1 6 , 0 . 0 , 0 . 0 , QD  ) 

CALL  LABLQ(ALPD,-20,QD,CRAPH,20) 

GO  TO  145 

143  CALL  XFSTQ(DISSR, 0.95,0. 18, 0.20, 0.0,0. O.QD)  I 

CALL  LABLQ(  ' SOUND  RAY  P ATH S '  , - 1 5 , Q D , C R A P H ,  1 5  ) 

145  C  LL  XFSTQ(DISSR, 0. 6 ,0. 18, 0.20, 0.0,0. O.QD) 
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IF(NAUT.EQ.O)  GO  TO  149 

CALL  LABLQ( ' RANGE  -  NAUTICAL  MILES ' , -2 2 , QD , GRAPH , 2 2 ) 

GO  TO  130 

149  CALL  LABLQ( ' RANGE  -  METERS -1 4 , QD , GRAPH , 1 4 ) 

150  CALL  XFSTQ(0 .0,0.0, 1 .0, 1 .0,0. 0,0.0 ,QD) 

CALL  AXISQ( RANL , NRAN ,0.0, RAN  INC , -0 . 13,-1 ,QD, GRAPH) 

CALL  DISPQ(GRAPH, 200.0) 

CALL  REMVQ ( GRAPH ) 

DELS V  -  FLOAT (NSV)*SVINC/SVL 
DELDEP  -  FLOAT (NDEP)*DEPINC/DEPL 
DELRAN  -  FLOAT(NRAN ) *RAN INC/ RANL 
C  DRAW  SOUND  VELOCITY  PROFILES 
DO  194  M  •  1 , NRSVPS 
AM  ■  M  -  1 

DISSV  -  RANL  +  SVL/2.0  +  1.0  +  AM*(SVL  +  1.0) 

CALL  XFSTQ( DISSV, 0.99 ,0. 18, 0.20,0. 0,0. 0,QD) 

IF(NRSVPS.LT.2)  GO  TO, 165 
CALL  EDINCH(M,TSVP(5) , 1  ) 

CALL  LABLQ(TSVP,-5,QD,GRAPH,5) 

GO  TO  166 

165  CALL  LABLQ( 'SVP'  ,-3 ,QD, GRAPH , 3  ) 

166  CALL  XFSTQ(DISSV, 0.55,0. 18, 0.20, 0.0, 0.0, QD) 

CALL  LABLQ( ' VELOCITY  -  M/ S EC ' , - 1 7 , QD , GRAPH , 1 7 ) 

DISVP  -  RANL  +  1.0  +  AM* ( S VL  +  1.0) 

CALL  XFSTQ( DISVP ,0.0,1.0,1.0,0.0,0.0,QD) 

CALL  AXISQ(SVL,NSV,SVMIN,SV INC, -0.1 3,-1 ,QD, GRAPH) 

DISSD  -  RANL  +  SVL  +  1.0  +  AM*(SVL  +  1.0) 

CALL  XFSTQ(DISSD,0.Q, 1.0,1.0,1.5*PI,0.0,QD) 

CALL  AXISQ(DEPL,NDEP,0.0,DEPINC,-0. 13,-1 , QD , GRAPH) 

DISDH  -  DISSD  +  0.52 
DISD  *  DEPL/2.0 

CALL  XFSTQ(DISDH,-DISD,0. 18,0.20,1.5*PI,0.0,QD) 

CALL  LABLQ( ' HEIGHT  -  METERS ',- 1 5 , QD , GRAPH , 1 5 ) 

DSVP  -  RANL  +  ( V ( M , 1 )  -  SVMIN)/DELSV  +  1.0  +  AM* ( SVL+1 . 0 ) 
CALL  XF STQ (DSVP, -0.02, 0.18,0. 20, 0.0,0.0,QD) 

CALL  LABLQ( '0'  ,-l ,QD, GRAPH, 1  ) 

CALL  ADPTQ(DSVP ,0 .0 , 1 .GRAPH) 

NSVP  -  NPSVP(M) 

DO  191  I  -  2, NSVP 

DSVP  -  RANL  +  ( V ( M , I )  -  SVMIN)/DELSV  +  1.0  +  AM*(SVL  +  1.0) 

DEPP  -  -  D ( M , I ) / DELDEP 

CALL  ADPTQ(DSVP , DEPP , 1 .GRAPH) 

IF(IGSVP(M,I) . LT  .  1 )  GO  TO  191 

CALL  XFSTQ(DSVP  ,DEPP  -  0 . 0 2 , 0 . 1 8 , 0 . 20 , 0 . 0 , 0 . 0 , QD ) 

CALL  LABLQ( '0'  ,-l ,QD,  GRAPH,  1  ) 

CALL  ADPTQ(DSVP , DEPP , 1 .GRAPH) 

191  CONTINUE 

CALL  DISPQ(GRAPH, 200.0) 

CALL  REMVQ(GRAPH) 

194  CONTINUE 
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CALL  ADPTQCRANL, 0.0,0, GRAPH) 

IF ( RANGE .LE.O.O;  GO  TO  19V 

CALL  ADPTQC RANGE/ DELRAN , -DEPL , 0 , GRAPH) 

CALL  ADPTQC RANGE /DELRAN ,0.0, I .GRAPH) 

199  CALL  XFSTQCO. 0,0.0, 1 .0, l .0, I . 5 *P I , 0 . 0 , QD ) 

CALL  AXISQCDEPL ,NDEP,0.0, DEPINC ,0. 13 ,-l ,QD .GRAPH) 
CALL  XFSTQ( -0 . 6  6 ,-DISD.O. 18,0.20, 1 .5*P 1,0.0, QD) 
CALL  LABLQ( 'HEIGHT  -  METERS 1 5 , QD , GRAPH , 1 5 ) 
IFCNRSVPS.LT. 2)  GO  TO  232 
READ (5,711)  (RSVP(I),  I  -  l.NRSVPS) 

DO  209  I  -  2 , NRSVPS 
RSV  -  RSVP ( I ) /DELRAN 

CALL  ADPTQCRSV,-  DEPL  -  1.0,0, GRAPH) 

CALL  ADPTQCRSV, 0.0,1, GRAPH) 

209  CONTINUE 

DISVP  -  DEPL  +  1.0 
NRSVPM  -  NRSVPS  -  1 
DO  217  I  -  1, NRSVPM 

RSV  -  0 .5*(RSVP( 1+1 )  +  RSVP( I) )/DELRAN 
CALL  EDINCH(I,TSVP(5),1) 

CALL  XFSTQCRSV , -DISVP, 0. 18, 0.20, 0.0, 0.0, QD) 

CALL  LABLQ(TSVP,-5 ,QD,GRAPH,5) 

217  CONTINUE 

CALL  EDINCH(NRSVPS ,TSVP(5)  ,  1  ) 

RSV  -  0 . 5 * ( RANL* DELRAN  +  RS VP ( NRSVPS ))/ DELRAN 
CALL  XFSTQCRSV, -DISVP, 0.1 6,0. 1 6 , 0 . 0 , 0 . 0 , QD ) 

CALL  LABLQ(TSVP,-5,QD .GRAPH, 5 ) 

CALL  DISPQCGRAPH, 200.0) 

CALL  REMVQ(GRAPH) 

IF(NAUT.EQ.O)  GO  TO  227 
WRITE(6,7l2)  TITLE 
GO  TO  228 

227  WRITE(6,713)  TITLE 

228  WRITE (6,714)  (I,RSVP(I),  I  -  l.NRSVPS) 
IF(NAUT.EQ.O)  GO  TO  232 

DO  231  I  -  l.NRSVPS 

231  RSVP ( I  )  -  6080.0*RSVP(I)/3.0 

232  IF ( NRBOT . LE . 0  )  GO  TO  264 

READ (5,715)  ( BDEP( I) , BRANC I)  ,  I  «  1, NRBOT) 

DO  911  1-1 , NRBOT 

911  BDEP ( I )-HTMAX-(BDEP( I )-HTMIN) 

IF ( METER . EQ  .  1 )  GO  TO  238 
LINE  ADDED  BDEP  STAYS  AS  METERS 
DO  237  I  -  1, NRBOT 
ID  -  328.084*BD£P(I) 

237  BDEP ( I  )  -  DFLOATC ID)/ 100 .0 

238  XD  »  200.0*BRAN( 1 )/DELRAN 
YD  -  -  200.0*BDEP(1 )/DELDEP 
CALL  DRAWQ ( XD , YD , 0 ) 

DO  245  I  -  2, NRBOT 
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XD  -  200.0*BRAN(I) / DEL RAN 
YD  -  -  200.0*BDEP(I)/DELDEP 
CALL  DRAWQ( XD ,  YD  ,  1  ) 

246  CONTINUE 

CALL  CDMPQ 

IF(NAUT.EQ.O)  GO  TO  250 
DO  24V  I  -  I , NRBOT 
24V  BRAN ( I )  -  60  80. 0* BRAN ( I ) / 3 . 0 
250  WRITE (6,716)  TITLE 
DEPMIN  «  BDEP(l) 

IM  -  NRBOT  -  1 
DO  260  I  -  1 , IM 
IP  -  I  +  1 

SLOPE(I)  -  (BDEP(I)  -  BDEP(IP))/( BRAN(  IP)  - 
XBRAN  -  3.0*BRAN(I)/6080.0 

WRITE (6, 71 7)  BDEP(I),XBRAN,BRAN(I),SLOPE(I) 
IF(DEPMIN.LT.BDEP(IP))  GO  TO  260 
DEPMIN  -  BDEP(IP) 

260  CONTINUE 

XBRAN  -  3 .0* BRAN (NRBOT)/ 6080 .0 

WRITE (6 ,717)  BDEP( NRBOT) .XBRAN, BRAN (NRBOT) 

GO  TO  265 

264  DEPMIN  -  100000.0 

265  IF(NAUT.EQ.O)  GO  TO  267 
DELR  -  6080.0*DELR/3.0 

267  NRSVPS  -  NRSVPS  +  1 

NPS VP ( NRS VP S )  -  NPSVP(NRS VPS-1 ) 

NRAY  -  1 

270  READ ( 5,718)  SD , NOSUR , NOBOT , RA , RMAX 
SD-HTMAX-(SD-HTMIN) 

IF(RMAX.LE.O .0)  GO  TO  611 
IF(NAUT.EQ.O)  GO  TO  274 
RMAX  -  6080 . 0* RMAX/ 3 . 0 
274  RSV P (NRSVPS )  -  RMAX  +  3000.0 
IF(METER.EQ.O)  GO  TO  278 
ID  -  328 . 084*SD 
C  SD  -  DFLOAT(ID)/100.0 

278  NSVP  -  NPSVP(l) 

DO  281  I  «  1 , NS  VP 
IF ( S D  -  D(1,I))286,284,281 
281  CONTINUE 

WRITE (6, 719)  SD,D(1,NSVP) 

GO  TO  270 
284  K  -  I 

GO  TO  311 
286  K  -  I 

J  -  I  -  1 
NR  -  1 

289  NPSVP(NR)  -  NPSVP(NR)  +  1 
NSVP  -  NPSVP(NR) 


BRAN ( I ) ) 
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GC  -  G (NR , J ) 

X  -  (SD  -  D  (  N  R  ,  J  ) )/ (  D  (  N  R  , I)  -  D(NR,J)) 

VC  -  V ( NR , J )  +  X*(V(NR, I)  -  V ( NR , J )  ) 

KB  -  I  +  1 
KOUNT  -  0 

DO  303  M  -  KB , NSVP 
MK  -  NSVP  -  KOUNT 
MN  -  MK  -  1 
D ( NR , MK )  -  D( NR, MN ) 

G ( NR , MK )  -  G( NR ,MN ) 

V(NR.MK)  -  V ( NR , MN ) 

IGSVP(NR.MK)  -  IGS VP ( NR , MN ) 

303  KOUNT  -  KOUNT  +  1 
D ( NR , I  )  -  SD 
G ( NR , I  )  -  GC 
V  (  NR  ,  I  )  =■  VC 
IGS VP ( NR , I )  -  0 
NR  -  NR  +  1 

IF(NR.GT.NRSVPS)  GO  TO  311 
IF(NPSVP(NR) .GE.K)  GO  TO  289 
311  KREST  *  0 
NOANG  -  0 
NOFIN  -  0 
NOM  -  0 
NR  -  1 

ANGLE ( 1 )  -  RA 
DEP(l)  -  SD 
R ( 1 )  -  0.0 
L  *  2 

NSVP  -  NPSVP(NR) 

THONE  -  0. 01745329252* RA 
CSTH  -  DCOS(THONE) 

SITH  -  DS IN ( THONE  ) 

TANTH  -  SITH/ CSTH 
VS  -  V( 1 ,K) 

SR  -  CSTH/VS 
Cl  -  G  C 1 , K- 1 ) 

G  1 2  -  G  (  1  ,  K  ) 

GO  TO  335 

330  IF(K.E;.  1)  GO  TO  333 

G I  -  FUNC(G(NR,K-1 ) ,R(L-1 ) ,RSVP(NR) ,RSVP(NR+i ) ,C(NR+1  ,K-1  )  ) 
IF(DABS(CI) .LT. 0.0002)  GI  -  0.0 
333  G I  2  -  FUNC(C(NR,K) ,R(L-1 ) ,RSVP(NR) ,RSVP( NR+1 ) , G ( NR+1 , K ) ) 
IF(DABS(C12) .LT. 0.0002)  GI2  -  0.0 

335  IF(SITH)348,336,338 

336  T  F  T.NPSVP(NR). AND. GI2.LT. 0.0)  GO  TO  348 
,E  .  1  .  OR .  C I  .  LE  .  0 . 0  )  GO  TO  366 

C  +  1 .0 
K  -  K  -  1 

IF(K.CT.O)  CO  TO  346 


338 
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IF(NOSUR.GT.O)  GO  TO  526 
K  -  K  +  1 
343  SITH  -  -  SITH 
TANTH  -  -  TANTH 
GO  TO  330 
346  GRAD  -  GI 
GO  TO  354 
348  C  -  -  1.0 
K  -  K  +  1 

IF(K.LE.NSVP)  GO  TO  353 
K  -  K  -  1 

IF(NOBOT)343,343,526 

353  GRAD  -  GI2 

354  IF(GRAD.NE.O.O)  GO  TO  373 
TATH  -  DABS (SITH/CSTH) 

IF(TATH.LT. 0.001)  GO  TO  366 
DD  -  DABS(D(NR,1C)  -  DEP(L-l)) 

DRFT  -  DD/TATH 

DR  -  DRFT 

SITH2  -  SITH 

CSTH2  -  CSTH 

TANTH2  -  TANTH 

ANGLE ( L )  -  ANGLE(L-l) 

DEP(L)  -  D ( NR , K ) 

GO  TO  401 

366  R(L)  -  R(L>1)  +  DELR 
DEP(L)  -  DEP(L-l) 

ANGLE ( L )  -  0.0 
SITH2  *  SITH 
CSTH2  -  CSTH 
TANTH2  -  TANTH 

IF ( DEP ( L )  -  DEPMIN)514,414,414 
373  VS  »  VS  +  GRAD* ( D ( NR , K)  -  DEP(L-l)) 
CSTH2  -  SR*V S 
B  -  1.0  -  CSTH2**2 
Y  -  DSQRT(DABS(B) ) 

IF ( Y . GE .0.001  )  GO  TO  382 

CSTH2  -  1.0 

SITH2  -  0.0 

TANTH2  -  0.0 

GO  TO  384 

382  IF(B.LT.O.O)  GO  TO  393 
SITH2  -  C* Y 

384  DR  -  DABS ( ( S ITH  -  SITH2)/ ( SR*CRAD ) ) 

ANGLE ( L )  -  57 .29577951* DA RSI N(SITH2) 
IF(DABS(ANGLE(L)  )  . LE  .  8  5 .0)  GO  TO  389 
NOANG  -  1 
GO  TO  527 

389  DEP ( L )  -  D ( NR , K ) 

R(L)  -  R(L-l)  +  DR 
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TANTH2  -  SITH2/CSTH2 
GO  TO  402 

393  KREST  -  1 

DDEPTH  -  (1.0  -  CSTH)/(SR*GRAD) 

VS  -  VS  -  GRAD* (  D ( NR , K )  -  DEP(L-l)) 

IC  -  C 

K  -  K  +  IC 

ANGLE (L)  -  0.0 

DEP(L)  -  D(NR, K)  +  DDEPTH 

DR  -  DABS((SITH/ (SR* GRAD))) 

401  R(L)  -  R(L-l)  +  DR 

402  IF(DEP(L) .GE.DEPMIN)  GO  TO  414 
L  -  L  +  l 

IF(KREST.EQ.O)  GO  TO  515 

405  DEP(L)  -  D(NR,K) 

ANGLE ( L )  -  -  ANGLE ( L-2 ) 

R ( L )  -  R(L-l)  +  DR 

SITH2  -  -  SITH 

CSTH2  -  CSTH 

TANTH2  -  -  TANTH 

IF (KREST . EQ . 1 )  KREST  -  0 

IF(DEP(L) .GE.DEPMIN)  GO  TO  414 

GO  TO  514 

414  DO  418  I  -  2.NRB0T 

IF(R(L-1 ) .GT.BRAN(I) )  GO  TO  418 
I BOT  »  I  -  1 
GO  TO  420 

418  CONTINUE 
GO  TO  514 

420  RC  -  DABS(BRAN(IBOT)  -  R(L) ) 

RC2  -  DABS(BRAN(IBOT+l )  -  R(L)  ) 

P  IF(RC.LT.0.1 .AND.DEP(L) . EQ . BDEP ( IBOT ) )  GO  TO  425 
IF(RC2  .GE .0. 1  )  GO  TO  427 
IF(DEP(L) .NE. BDEP(IB0T+1 ) )  GO  TO  427 

425  NOM  -  1 
GO  TO  527 

427  ANGX  -  0.01745329252* ANGLE(L-1  ) 

SITHX  -  DSIN(ANGX) 

CSTHX  -  DCOS(ANGX) 

TANTHX  -  SITHX/CSTHX 
IF(SITHX.LT.O.O)  GO  TO  435 
IF(DEP(L~1 ) .GT .BDEP(IBOT)  )  GO  TO  434 
IF(DEP(L-1).LT.BDEP(IB0T+1))  GO  TO  511 

434  IF(SLOPE( IBOT) ) 5 1 1 , 51 1 ,441 

435  IF(DEP(L) . GT . BDEP ( IBOT)  )  GO  TO  441 
IF(DEP(L) .GT.BDEP( IBOT+1 ) )  GO  TO  441 

IF(SLOPE(IBOT).EQ.O.O.AND.DEP(L).EQ.BDEP(IBOT))  GO  TO  477 
IF(ANGLE(L)  .NE.O.O)  GO  TO  511 
KREST  -  1 
GO  TO  511 
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441  DM  -  BDEP(IBOT)  -  DEP(L-l)  -  SLOPE ( IBOT )*( R ( L-l )  -  BRAN ( IBOT ) 
IF(GRAD.NE.O.O)  GO  TO  446 
IF(SLOPE(IBOT) .EQ.TANTHX)  GO  TO  511 
DELX  -  -  DM/ (TANTHX  -  SLOPE(IBOT)) 

GO  TO  458 

446  GAM  -  1 .0/(SR*GRAD) 

ALPH  -  -  GAM*S ITHX 
BETA  -  -  GAM*CSTHX 
AQ  -  SLOPE( IBOT)**2  +  1.0 

BQ  -  2 . 0* ( -  SLOPE(IBOT)*(DM  -  BETA)  -  ALPH) 

CQ  -  ALPH**2  +  (DM  -  BETA)**2  -  GAM**2 

DISC  -  BQ**2  -  4.0*AQ*CQ 

IF (DISC. LT. 0.0)  GO  TO  511 

IF ( GRAD .GT.0.0)  GO  TO  457 

DELX  -  (-  BQ  +  DSQRT(DISC))/(2.0*AQ) 

GO  TO  458 

457  DELX  -  (-  BQ  -  DSQRT(DISC))/(2.0*AQ) 

458  IF ( DABS ( DELX  -  DR).LE.0.1)  GO  TO  474 
IF(DELX.LT.O.O)  GO  TO  511 

XP  -  R(L-l)  +  DELX 

IF(XP.LT. BRAN ( IBOT )  )  GO  TO  511 

IF(XP.GT.BRAN( IBOT+1 )  )  GO  TO  511 

IF ( XP . GT . R ( L  )  )  GO  TO  511 

DELY  -  -  DELX*SLOPE( IBOT)  +  DM 

R(L)  -  XP 

DEP(L)  -  DEP(L-l)  +  DELY 

VS  -  VS  -  GRAD* (D(NR ,K)  -  DEP(L) ) 

CSTH2  -  SR* VS 

IF(CSTH2.GT  .  1 .0)  GO  TO  4/5 
IF(ANGLE( L-l ) .EQ.O .0 )  C  -  -  C 
SITH2  -  C*DSQRT( 1.0  -  CSTH2**2) 

TANTH2  -  SITH2/CSTH2 

ANGLE ( L )  -  57 .29577951* DA RSIN(SITH2) 

474  IF(DABS(ANGLE(L) ) .LE.85.0)  GO  TO  477 

475  NOANG  -  1 
GO  TO  527 

*77  IF(R(L) .GT.RMAX)  GO  TO  527 
IF( NOBOT.EQ. 1 )  GO  TO  527 
L  -  L  +  1 
DEP(L)  -  DEP(L-l) 

R(L)  -  R(L-l) 

ANGLE ( L  )  -  114. 59156* DA TAN(SLOPE(IBOT))  -  ANGLE ( L-l ) 
IF(DABS(ANGLE(L)  )  . LE.85.0)  GO  TO  486 
NOANG  -  1 
GO  TO  527 

486  THONE  -  0.01745329* AN GLE(L) 

SITH  -  DS IN ( THONE  ) 

CSTH  -  DCOS ( THONE ) 

TANTH  -  S ITH/ CSTH 

IF ( GRAD .EQ.0.0)  GO  TO  497 
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IF(KREST.EQ.O)  GO  TO  *95 
KREST  *  0 

IF(SITH2.LE.O.O.AND.S1TH.GT.O.O)  K  -  K  +  1 
GO  TO  497 

495  IF(SITH.GT.0.0 .AND .SITH2.GT.0.0)  K  -  K  +  1 
IF(SITH.LT. 0.0. AND. SITH2.LT. 0.0)  K  -  K  -  1 
497  L  «  L  +  1 

SR  -  CSTH/VS 

IF ( CRAD .NE.O.O)  GO  TO  330 

IF(SITH.GT. 0.0. AND.SITH2.LT. 0.0)  K  »  K  -  1 
DD  -  DEP(L-l)  -  D (NR, K) 

DRFT  -  DABS(DD/TANTH) 

DR  -  DRFT 
DEP(L)  ■  D( NR , K ) 

R(L)  -  R(L-I)  +  DR 
ANGLE ( L  )  -  ANGLE(L-I) 

SITH2  -  SITIt 
CSTH2  «  CSTH 
TANTH2  -  TANTH 
GO  TO  514 

511  IF(R(L).LT.BRAN(IB0T+1))  GO  TO  514 
I BOT  -  I BOT  +  1 
GO  TO  42/ 

514  L  -  L  +  1 

515  IF ( L . LT .4000 )  GO  TO  518 
NOFIN  -  l 

GO  TO  52b 

518  IF(R(L-1 ) .GE.RMAX)  GO  TO  526 
IF ( KREST . EQ .  1 )  GO  TO  405 
SITH  -  SITH2 
CSTH  -  CSTH2 
TANTH  -  TANTH2 

IF(R(L-1 ) .LE.RSVPCNR+l ) )  GO  TO  330 
NR  -  NR  +  1 

IF(NR.LT.NRSVPS )  GO  TO  330 

526  L  -  L  -  1 

527  IF(NAUT.EQ.O)  GO  TO  530 
DO  529  I  -  l,L 

529  R(I)  -  3  »  0*R( I )/6080 .0 

530  IF(NOPR.GT.O)  GO  TO  549 
IB  -  1 

NPAGE  •  1 

533  WRITE( 6,720)  NPAGE , TITLE , SD , RA 
IF(NAUT.EQ.O)  GO  TO  537 
WRITE(6,721 ) 

GO  TO  538 

537  WRIT£(6 , 722  ) 

538  IL  -  IB  +  43 
IF(IL.LE.L)  GO  TO  541 
IL  -  L 


I 


4 


4 


\ 


4 
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541  WRITE (6 ,  723 )  (DE P ( I ) , R ( I ) , ANGLE ( I )  ,  I  -  IB,IL) 

IB  -  IL  +  1 
IF ( IB  -  L ) 544 , 547,552 
544  WRITE (6 , 710) 

NPAGE  -  NPACE  +  1 
GO  TO  533 

547  WRITE (6,723)  DEP (L ) ,R (L ) .ANGLE (L ) 

GO  TO  552 

549  IF (NOMO.GT .  1  )  GO  TO  551 
WRITE (6 , 724 )  TITLE, SD 

551  NOMO  -  NOMO  +  1 

552  IF(NOFIN.EQ. 0)  GO  TO  554 
WR ITE(6, 725)  RA 

554  IF(NOM.EQ.O)  GO  TO  556 
WRITE ( 6 , 726)  RA 
556  IF(NOANG.EQ.O)  GO  TO  558 
WRITE(6, 727)  RA 
558  WRITE (6 , 728)  RA , L 

RMAXL  -  RMAX  +  0.5*RANINC 
IF ( R (L ) . LE . RMAXL )  GO  TO  566 

561  DEP(L)  -  DEP(L-l)  +  (DEP(L)  -  DE P (L-l ) ) * (RMAX  -  R(L-1))/(R(L) 
1  R (L-l ) ) 

R (L )  -  RMAX 

IF(R(L-1 ) .LE.RMAX)  GO  TO  566 
L  -  L  -  1 
GO  TO  561 

566  IF (NRAY . LT . 2  )  GO  TO  571 
IF (NOBOT . LT.  1  )  GO  TO  569 
L  -  L  -  1 
569  NRAY  -  1 
GO  TO  572 

571  NRAY  -  2 

572  PRAN  -  R (L ) /DELRAN  +  0.05 
PDEP  -  -  DE  P (L ) /DELDE  P  -  0.05 
IF (NRAY . LT . 2 )  GO  TO  601 

57£  XD  -  200. 0*R ( 1 ) /DELRAN 

YD  -  -  200. 0*DEP (  1  ) /DELDEP 
CALL  DRAWQ ( XD , YD , 0 ) 

IB  -  2 

IF(L.CT. 990)  GO  TO  582 
IL  -  L 
GO  TO  583 

582  IL  -  990 

583  DO  587  I  -  IB.IL 

XD  -  200. 0*R(I ) /DELRAN 
YD  -  -  200. 0*DEP(I )/DELDEP 
CALL  DRAWQ (XD, YD, 1) 

587  CONTINUE 
CALL  GDMPQ 
IB  -  IL  +  l 
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IF(IB.GT.L)  GO  TO  595 
IL  -  IL  +  990 
IF(IL.LE.L)  GO  TO  583 
IL  -  L 
GO  TO  583 

595  IF (NRAY. LT . 2  )  GO  TO  270 

596  CALL  XFSTQ(PRAN, PDEP, 0. 1 l, 0. 13, 0.0,0. 0,QD) 

CALL  NMBRQ(RA, I,  1 , QD , GRAPH) 

CALL  DISPQ(GRAPH, 200. 0) 

CALL  REM VQ (GRAPH ) 

GO  TO  (575, 270) , NRAY 
601  N  -  L/2 

DO  609  I  -  1 , N 
K  -  L  -  I  +  l 
TR  -  R (  I  ) 

TD  -  DEP(I) 

MI)  -  R(K) 

DEP(I)  -  DEP(K) 

MK)  -  TR 
609  DEP(K)  -  TD 
GO  TO  596 

611  READ (5 , 704  )  DIST 

IF(DIST.LE.O.O)  GO  TO  615 
CALL  ADPTQCDIST, 0. 0, 0, GRAPH) 

GO  TO  616 

615  CALL  ADPTQ(0. 0,0. 0,0, GRAPH) 

616  CALL  DISPQ(GRAPH, 200.0) 

CALL  MOO  VQ ( 200 . 0*DIST , 0. 0) 

CALL  REM  VQ (GRAPH ) 

NOMO  -  1 

IF(DIST.GT.O.O)  GO  TO  10 
WRITE ( 6 , 729) 

STOP 

701  FORMAT ( 4  OA  1 ) 

702  FORMAT (I  1, 13, 21 1 , F4 . 0 , I  2 , 311 , 3 (I  2 , F8 . 2 , F5 .  1  )  , 2F  1  0.  4  ) 

703  FORMAT (711) 

704  FORMAT (4F10.4,2I1) 

705  FORMAT ( IH l , 1 2 1 X , ' PAGE ' , 1 3 , / IHO, 88X , 40A 1 , / 8  9X'NUMB  ER  OF  POINTS  IN', 
1  '  SVP',12,'  -  ',14) 

706  FORMAT(89X, 'VELOCITIES  COMPUTED') 

707  FORMAT ( 1H0, 61X, 'SVP  ' , 1 1 , / 1H0, 33X, ' DEPTH  TEMPERATURE ', 3X, 

1  'SALINITY  VELOCITY  VELOCITY  GRADIENT' /35X, ' (FT) ', 6X, 

2  '(DEG  F)' ,6X,' (?PT)' , 5X , ' (FT/SEC)' ,6X,' (FT/SEC/FT)'  /) 

708  FORMAT ( IH  , 2 9X , F9 .  1 , 2F 1 2 . 2 , F 1 3 . 3 , / 8 OX , F  1  2 . 7 ) 

709  FORMAT ( IH  , 2 9X , F9 .  1 , 2 5X , F 1 2  .  3 , / 80X , F 1 2 .  7 ) 

710  FORMAT (l HO, 58X,' (CONTINUED)' ) 

711  FORMAT (6F10.4) 

712  FORMAT ( IH  l, 88X,  40A 1, /1H0,58X,'SVP  RANGES' /1H0,55X,'SV?',8X, 

1  'RANGE'/68X,  '  (NM) ' ) 

713  FORMAT (1H1,88X,40A 1 , / IH 0 , 58X , ' S VP  RANGES '/1H0, 55X, 'SVP' , 8X, 
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1  'RANGE'/67X,' (YDS)') 

714  FORMAT (1H0, 55X, 12, F  14  .  1) 

715  FORMAT(2F 10. 4) 

716  FORMAT(1H1,88X,40A1,/1HO,47X, 'BOTTOM  DEPTHS,  RANGES  AND  SLOPES' 

1  /1H0, 44X, 'DEPTH' , 1 2X ,' RANGE ', 1 2X ,' S LOPE ' /4 6X, '( FT )', 7X ,' (NM )', 7X , 

2  '(YDS)'  /) 

717  FORMAT (40X,F11.  l,F10.l,F12.1,/7  3X,F11.4) 

718  F ORM AT (F8.2,2I1,2F10.4) 

719  FORMAT( 1H 1 , 30X, 'REQUESTED  DEPTH  0F',F7.1,'  FT.  IS  GREATER  THAN  ', 

1  'LAST  GIVEN  SVP  DEPTH  0F',F7.1,'  FT.') 

720  F0RMAT(IH1, 121X, 'PAGE', 13, / 1H0 , 88X , 40A 1 , / 8 9X , 'SOURCE  DEPTH  -', 

1  F 8 . 1 , '  FT. '/89X, 'INITIAL  ANGLE  -',F8.3,'  DEG. ' /1H0, 56X, ' SOUND  RAY 

2  PATH'/lHO, 4 9X, 'DEPTH' , 7X, 'RANGE ' , 7X, 'ANGLE ' ) 

721  F0RMAT(51X, '  (FT) '  , 8X,  '  (NM) ' , 7X, '  (DEG) '  /) 

722  FORMAT (5 IX, ' (FT) ' , 7X, ' (YDS ) ' , 7X, ' (DEG) '  /) 

723  FORMAT(45X,F10. 1,F12.1,F12.2) 

724  FORMAT (1H1,88X,40A1,/89X, 'SOURCE  DEPTH  «',F8.1,'  FT. ' /1H0, 49X, 

1  'NUMBER  OF  POINTS  IN  RAY  PATHS  AND ' / 4  OX , ' RAYS  THAT  TERMINATE  ', 

2  'BEFORE  REACHING  DESIRED  RANGE') 

725  FORMAT (1H0,40X,F5. 1,'  DEG.  RAY  TERMINATED.  MORE  THAN  4000  POINTS') 

726  FORMAT (1H0,40X,F5. 1,'  DEG.  RAY  TERMINATED.  HIT  BOTTOM  BREAK.') 

727  FORMAT ( IHO , 40X , F5 • l,'  DEG.  RAY  TERMINATED.  ANGLE  GREATER  THAN  85', 
1  '  DEG.  '  ) 

728  FORMAT(1HO, 40X, 'NUMBER  OF  POINTS  IN',F5.1,'  DEG.  RAY  -',15) 

729  FORMAT ( 1H 1 , 5 OX , ' RUN  COMPLETED') 

END 


Appendix  C 


The  input  for  the  eigenray  and  ray-tracing  programs  is  described 
here.  First,  the  input  to  the  eigenray  routine  contained  in  Appendix 
A  is  given. 


Eigenray  Program  Operation  -  Sequence  of  Data  Cards 
Columns 

1.  Title  card 

1-40  title 

2.  Sound  velocity  profile  control  card 

0  meters 

1  NDRT  -  depth  given  in 

1  feet 

0  m/sec 

2  NVFT  *  sound  velocity  in 

1  ft/sec 


Format 

10A4 

911 


NTF 


temperature  given  in  degrees  ^ 


0  m/sec 

NWV  *  wind  velocity  given  in 

1  ft/sec 

0  meters 

NRFT  *  range  given  in 

i  feet 

0  calculated 

IVEL  *  sound  velocity 

1  given 

0  not  given 

ITEMP  »  temperature 

1  given 

0  not  given 

IWV  -  wind  velocity 

1  given 

0  not  given 

IRHM  *  relative  humidity 

1  given 


3.  Sound  velocity  profile  data  cards  -  must  be  given  as  4F10.4,lx,Il 
depths  from  highest  level  to  ground  surface 

I- 10  DEP  -  depth 

II- 20  TEMP  -  temperature 

21-30  VEL  -  sound  velocity 

31-40  WV  -  wind  velocity 

42  NOMO  *1  to  indicate  last  SVP  data  card 
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* 


4. 

Ground  loss  card 

11 ,4X,2F5*1 

1 

NBL  -  number  of  runs  with  different 

ground  losses 

6-10 

POR  -  porocity  of  ground 

11-15 

BL  -  ground  loss  coefficient 

5 

Time  integration 

wii.dow  card 

F10.3 

1-10 

TWIN  -  time  window 

6. 

Ray  path  parameter  card 

7F10.3 

1-10 

SD  -  source  depth 

11-20 

TD  -  target  depth 

21-30 

RANGE  -  range  from  source  to  receiver 

31-40 

ANGMAX  -  maximum  initial  ray  angle 

in  degrees 

41-50 

ANGMIN  -  minimum  initial  ray  angle 

51-60 

FREQ  -  frequency 

61-70 

RHM  -  relative  humidity 

All  input  must  be  consistant  with  the  units  specified  in  the 
control  card  (2). 

The  input  for  the  ray-tracing  routine  is  given  here.  It  is  noted 
that  the  ray-tracing  routine  graphics  are  system  dependent.  T’.e  data 
to  be  plotted  is  output  using  the  following  packages:  XFSTQ,LABLQ, 
AXISQ,  DISPQ,  REMVQ,  ED INCH,  ADPTQ,  DRAWQ,  GDMPQ ,  NMBRQ,  MOOVQ. 

These  are  packages  available  in  the  PSU  computer  center's  accessible 
library.  The  output  is  then  plotted  onto  a  Tektronix  4662  plotter 
using  the  package  CONTK. 

Ray-Tracing  Program  Operation  -  Sequence  of  Data  Cards 


Columns 

Format 

Title  card 

1-40 

title 

40A1 

AD-A108  626  PENNSYLVANIA  STATE  UNIV  UNIVERSITY  PARK  NOISE  CONTROL  LAB  p/5  ?0/l 

RAY  TRACING  TECHNIQUES  -  DERIVATION  ANO  APPLICATION  TO  ATMOSPHE-- ETC (U) 
JAN  60  S  0  ROTH  DAAK07-60-C-0001 

NL 


UNCLASSIFIED 
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t 

i 


f 


2*  Sound  velocity  profile  control  card 

1  NRSVP  -  number  of  SVP's  -  maximum  II 

of  5 

2-4  NRBOT  -  number  of  ground  height  13 

coordinates  -  maximum  of  130 
minimum  of  2 

5  METER  *  1  all  heights  given  in  meters  II 

6  NAUT  -  blank  all  ranges  given  in  meters  II 

7-10  DELR  -  distance  ray  travels  when  ray  F4.0 

angle  is  zero  and  there  is  no 
refraction 

11-12  NOUT  -  7)  12 

blank  -  uses  SVP's  from  previous  run 
13  NEWS  VP  -  II 

1  -  new  SVP's 


blank  prints 

14  NOPR  •  ray  paths  II 

1  does  not  print 

blank  "SOUND  RAY  PATHS"  II 

15  NOAD  *  writes  on  plots 

1  "SOUND  RAY  PATHS  AND 

INTENSITY  CONTOURS" 
and  Alpha  headings 

16-17  NRAN  -  number  of  divisions  on  range  scale  12 

18-25  RANINC  -  Increment  on  range  scale  -  meters  F8.2 

26-30  RANL  -  length  in  inches  of  range  scale  F5.1 

31-32  NDEP  -  number  of  divisions  on  height  scale  12 

33-40  DEPINC  -  increment  on  height  scale  -  meters  F8.2 

41-45  DEPL  -  length  in  inches  of  height  scale  F5.1 

46-47  NSV  -  number  of  divisions  on  SVP  scale  12 

48-55  SVINC  -  increment  on  SVP  scale  -  m/sec  F8.2 

56-60  SVL  -  length  in  Inches  of  SVP  scale  F5.2 

61-70  SVMIN  -  minimum  value  on  SVP  scale  F10.4 

(generally  335  m/sec) 

71-80  RANGE  -  range  to  reference  line  in  meters  F10.4 

(if  blank-no  reference  line  printed) 

3.  Alpha  heading  card  (omitted  if  NOAD  Is  blank) 

1-20  ALPD  -  graph  heading  20A1 


V 


■  +  -4  '• 


*  *  - 
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One  sec  of  cards  4  and  5  must  be  given  for  each  sound  velocity  profile 

4,  SVP  control  card  711 

1  NODEP  *  blank 

2  NODFT  -  1 

blank  temperature  given 

3  NOTEMP  - 

1  sound  velocity  in  still  air  given 

4  NOTF  -  1 

5  NOVEL  -  blank 

6  NOVFT  -  1 

7  NO SAL  -  blank 

5.  Sound  velocity  profile  cards  -  maximum  of  250  4F10.4*2I1 

for  each  SVP 

1-10  D  -  height  meters 

11-20  TF  -  temperature  Celsius 

21-30  WV  -  wind  velocity  m/sec 

31-40  SAL  -  blank 

blank  interpolated 

41  IGSVP  -  SVP  value 

1  given 

blank  not  last  SVP  card 

42  NOMO  -  indicates 

1  last  SVP  card 

6*  Location  of  SVP's  card  (first  always  0.0)  5F10.4 

(omitted  if  only  one  SVP) 

1-10  RSVP(l)  -  0.0 

11-20  RSVP(N)  -  locates  SVP  meters 

etcetera 

/.  Ground  height  coordinate  cards  (must  be  NRBOT  cards*  2F10.4 

value  given  on  card  2)  (must  be  In  order  of 
increasing  range) 

1-10  B0EP  -  height  meters 

11-20  BRAN  -  range  meters 
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8.  Ray  parameter  card*  F8«2,2Il,2F10.4 

1-8  SO  -  aourca  height  meter* 

9  NOSUR  -  1 

10  NOBOT  »  blank 

11-20  RA  -  Initial  angle  in  degree* 

21*30  RMAX  *  maximum  range  in  meter* 


9 .  Blank  card 


DATE 

FILMED 


DTIC 


